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

TL;DR

This study uses molecular dynamics to explore ferroelectric glass behavior in spheroidal dipolar systems with impurities, revealing polar nanoregions, heterogeneity, and nonergodic responses similar to relaxor ferroelectrics.

Contribution

It introduces a detailed simulation of dipolar glasses with spheroidal particles and impurities, elucidating the formation of polar nanoregions and nonergodic phenomena.

Findings

Formation of polar nanoregions (PNRs) at low temperatures.

Heterogeneous impurity clustering affects dielectric response.

Observation of nonergodic polarization and strain behavior.

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 KCN$_x$KBr$_{1-x}$), 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 $T$ depending on the size of the dipole moment. They then form frozen PNRs as $T$ 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 (ZFC$/$FC), 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.

Equations

$$N_1=V{n}_1=(1-c)N,N_2=V{n}_2=cN,$$

$$U=U_{\mathrm{LJ}}+U_{\mathrm{w}}+U_{\mathrm{d}}.$$

$$U_{\rm LJ}=2\epsilon\sum_{i\neq j}\bigg{[}(1+A_{ij})\frac{\sigma_{\alpha\beta}^{12}}{r_{ij}^{12}}-\frac{\sigma_{\alpha\beta}^{6}}{r_{ij}^{6}}\bigg{]}.$$

$$A_{ij}=\delta_{\alpha 1}{\eta}({\mbox{\boldmath$n$}}_{i}\cdot{\mbox{\boldmath$r$}}_{ij}/r_{ij})^{2}+\delta_{\beta 1}{\eta}({\mbox{\boldmath$n$}}_{j}\cdot{\mbox{\boldmath$r$}}_{ij}/r_{ij})^{2},$$

$${\sigma }_2/{\sigma }_1=1.1,\eta =\mathrm{1.2.}$$

$$U_{\mathrm{w}}=w\sum _i[\exp (-z_i/{\xi }_{\mathrm{w}})+\exp (-(H-z_i)/{\xi }_{\mathrm{w}})],$$

$${\mbox{\boldmath$\mu$}}_{i}=(\mu_{xi},\mu_{yi},\mu_{zi})=\mu_{0}{\mbox{\boldmath$n$}}_{i},$$

$$\Phi (x,y,0)=0,\Phi (x,y,H)=-\Delta \Phi =-E_{\mathrm{a}}H,$$

$$\bar{\mbox{\boldmath$\mu$}}_{i}=(-\mu_{xi},-\mu_{yi},\mu_{zi}),$$

$$\displaystyle\Phi({\mbox{\boldmath$r$}})$$

$$\displaystyle U_{\rm d}=\frac{1}{2}\sum_{{\mbox{\boldmath$h$}}}{\sum_{i\in 1,j\in 1}}^{\prime}{\mbox{\boldmath$\mu$}}_{i}\cdot{\stackrel{{\scriptstyle\leftrightarrow}}{{\mbox{\boldmath$\cal T$}}}}({\mbox{\boldmath$r$}}_{ij}+{\mbox{\boldmath$h$}})\cdot{\mbox{\boldmath$\mu$}}_{j}$$

$$\displaystyle+\frac{1}{2}\sum_{{\mbox{\boldmath$h$}}}{\sum_{i\in 1,j\in 1}}{\mbox{\boldmath$\mu$}}_{i}\cdot{\stackrel{{\scriptstyle\leftrightarrow}}{{\mbox{\boldmath$\cal T$}}}}({\bar{\mbox{\boldmath$r$}}}_{ij}+{\mbox{\boldmath$h$}})\cdot{\bar{\mbox{\boldmath$\mu$}}}_{j}-E_{\rm a}M_{z}.$$

$${\mbox{\boldmath$M$}}=(M_{x},M_{y},M_{z})=\sum_{i}{\mbox{\boldmath$\mu$}}_{i}.$$

$${\mbox{\boldmath$E$}}_{i}=-\partial U_{\rm d}/\partial{\mbox{\boldmath$\mu$}}_{i}.$$

$${\mbox{\boldmath$E$}}_{i}=-\lim_{{\mbox{\boldmath$r$}}\to{\mbox{\boldmath$r$}}_{i}}\nabla[\Phi({\mbox{\boldmath$r$}})-{\mbox{\boldmath$g$}}({\mbox{\boldmath$r$}}-{\mbox{\boldmath$r$}}_{i})\cdot{\mbox{\boldmath$\mu$}}_{i}].$$

$$\mathcal{H}={\mathcal{H}}_0-M_z{E}_{\mathrm{a}},$$

$$\frac{\partial }{\partial E_{\mathrm{a}}}⟨\mathcal{A}{⟩}_{\mathrm{e}}=\frac{1}{k_{\mathrm{B}}T}⟨\mathcal{A}\delta M_z{⟩}_{\mathrm{e}},$$

$${\chi }_{\mathrm{dif}}=\frac{\partial P_z}{\partial E_{\mathrm{a}}}=\frac{1}{V{k}_{\mathrm{B}}T}⟨(\delta M_z{)}^2{⟩}_{\mathrm{e}}.$$

$${\cal K}=\frac{1}{2}\sum_{i}m|{\dot{\mbox{\boldmath$r$}}}_{i}|^{2}+\frac{1}{2}\sum_{i\in 1}I_{1}|{\dot{\mbox{\boldmath$n$}}}_{i}|^{2},$$

$$m{\ddot{\mbox{\boldmath$r$}}}_{i}=-{\partial U}/{\partial{\mbox{\boldmath$r$}}_{i}},$$

$$I_{1}({\ddot{\mbox{\boldmath$n$}}}_{i}+|{\dot{\mbox{\boldmath$n$}}}_{i}|^{2}{\mbox{\boldmath$n$}}_{i})=(\stackrel{{\scriptstyle\leftrightarrow}}{{\mbox{\boldmath$I$}}}-{\mbox{\boldmath$n$}}_{i}{\mbox{\boldmath$n$}}_{i})\cdot\mu_{0}{\mbox{\boldmath$E$}}_{i}^{\rm eff},$$

$${\mbox{\boldmath$E$}}_{i}^{\rm eff}={\mbox{\boldmath$E$}}_{i}+{\mbox{\boldmath$E$}}_{i}^{\rm ste},$$

$${\mbox{\boldmath$E$}}_{i}^{\rm ste}=-({8\varepsilon{\eta}}/{\mu_{0}})\sum_{j\neq i}({\sigma}_{\alpha\beta}^{12}/r_{ij}^{14})({\mbox{\boldmath$n$}}_{i}\cdot{{\mbox{\boldmath$r$}}}_{ij}){\mbox{\boldmath$r$}}_{ij}.$$

$$t_0={\sigma }_1(m/ϵ{)}^{1/2}.$$

$$Q_{\ell i}=\sum_{j\in{\rm neighbor}}P_{\ell}({\mbox{\boldmath$n$}}_{i}\cdot{\mbox{\boldmath$n$}}_{j})/Z_{i}\quad(\ell=1,2),$$

$$⟨(\delta Q_{\ell }{)}^2⟩=\sum _{i\in 1}(Q_{\ell i}(t_0)-⟨Q_{\ell }⟩{)}^2/N_1(\ell =1,2).$$

$$C_{1}(t)=\sum_{i\in 1}\langle{{\mbox{\boldmath$n$}}_{i}(t_{0})\cdot{\mbox{\boldmath$n$}}_{i}(t_{0}+t)}\rangle/N_{1},$$

$$C_1({\tau }_1)=0.1,$$

$${\tau }_1={\tau }_{10}\exp [D_1{T}_1/(T-T_1)].$$

$$\omega =T_{\mathrm{m}}^{-1}(T)={\omega }_{\mathrm{m}}(T),$$

$$V{k}_{\mathrm{B}}T{\chi }_{\mathrm{fl}}={⟨(\delta M_z{)}^2⟩}_{\mathrm{time}}={⟨M_z^2⟩}_{\mathrm{time}}-{⟨M_z⟩}_{\mathrm{time}}^2.$$

Full text

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 Gao ; Wei ; Weis ; Ta ; Tao ; Ayton ; Hent0 ; Hent ; Dij0 ; Groh . Such spherical dipoles form various noncubic crystals in ferroelectric phases Groh ; Dij0 . Ferroelectriciity was also studied in positionally disordered dipolar solidsAyton . Recently, Johnson et al.Johnson0 ; Johnson1 have investigated a ferroelectric transition of spheroidal particles with a dipole moment $\mu_{0}$ parallel to the spheroidal axis. They found that the static dielectric constant increases up to $10^{2}-10^{3}$ with increasing $\mu_{0}$ 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 disorderVug ; Bl ; Kob , which take place over a wide temperature range without long-range dipolar order. Notable examples are relaxors Smo ; Klee ; Klee1 ; Bokov3 ; Bokov ; Bl ; Cr ; Cowley ; Samara such as Pb(Mg1/3Nb2/3*)O3 (PMN), where the random distribution of Mg2+ and Nb5+ at B sites yields quenched random fields at Pb2+ sites Setter ; West . In relaxors, temperature-dependence of the optic index of refraction suggested appearance of mesoscopic polarization heterogeneities Burns , called polar nanoregions (PNRs). They are enhanced at relatively high $T$ as near-critical fluctuations and are frozen at lower $T$. It is widely believed that these PNRs give rise to a broad peak in the dielectric permittivity $\varepsilon^{\prime}$ as a function of $T$ Setter ; Smo ; Klee ; Klee1 ; Bl ; Cr ; Cr1 ; Bokov ; Bokov3 ; Cowley ; Samara ; Ra . They have been detected by neutron and x-ray scattering Wel ; Stock ; Eg ; Shi1 ; Shi2 ; Cowley and visualized by transmission electron microscopy Ban ; Mori ; Uesu and piezoresponse force microscopy Bokov2 ; Bokov4 ; Ann . Strong correlations have also been found between the PNRs and the compositional heterogeneity of the B site cations Wel ; Wel1 ; Ban ; Setter ; Bokov4 ; Mori ; Perrin ; Jin ; i3 ; i6 .

Relaxor behaviors also appear in other disordered dipolar systemsBl ; Kob ; Vug . In particular, orientational glass has long been studied in mixed crystals such as KCNxKBr1-x or KxLi1-xTaO3 Vug ; Tou ; Ho1 ; Ho2 ; Loidl3 ; qua ; ori ; Binder ; anti-Loidl ; Yokota ; Mertz ; F3 , 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 $T$, 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 heterogeneitiesqua ; Heuer ; Takae-ori ; EPL . Some polymers also undergo ferroelectric transitions due to alignment of permanent dipoles Bl ; Lov ; Furukawa ; Furu1 . In particular, poly(vinylidene fluoride-trifluoroethylene) copolymersZhang ; Z1 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 F5 , 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 $T$.

In their molecular dynamics simulation of relaxors, Burton et al.Burton1 ; Burton2 ; Burton3 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 apply ; Klapp , which has been used to study water between electrodesHautman ; Takae1 ; Takae2 . To prepare a mixed crystal, we cool a liquid mixture from high $T$; then, our impurity distribution at low $T$ is naturally formed during crystallization Takae-ori ; EPL .

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 ZFC$/$FC (zero-field-cooling and field-cooling) temperature changes. In the latter, nonergodicity of glass is demonstrated, so its experiments have been performed in spin glassMydosh ; F0 ; F1 , relaxorsWest ; Uesu ; F4 , orientational glassHo1 ; F3 , relaxor-like polymers Z1 , and strain glassF5 .

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 $z=0$ and $H$ to apply electric field to the dipoles. The periodic boundary condition is imposed along the $x$ and $y$ axes with period $L$. Thus, the particles are in a $L\times L\times H$ cell with volume $V=L^{2}H$.

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

[TABLE]

where the total particle number $N$ is set equal to $8000$. Their positions are written as {\mbox{\boldmathr}}_{i}=(x_{i},y_{i},z_{i}) ($1\leq i\leq N$). The long axes of the spheroidal particles are denoted by unit vectors {\mbox{\boldmathn}}_{i}=(n_{xi},n_{yi},n_{zi}) ($1\leq i\leq N_{1}$).

II.1 Potential energy

The total potential energy $U$ is expressed as

[TABLE]

Here, $U_{\rm LJ}$ is the sum of modified Lennard-Jones potentials between particles $i\in\alpha$ and $j\in\beta$ ($\alpha,\beta=1,2$),

[TABLE]

where r_{ij}=|{\mbox{\boldmathr}}_{i}-{\mbox{\boldmathr}}_{j}|, $\epsilon$ is the characteristic interparticle energy, and $\sigma_{\alpha\beta}=(\sigma_{\alpha}+\sigma_{\beta})/2$ in terms of the particle lengths $\sigma_{1}$ and $\sigma_{2}$. The factor $A_{ij}$ depends on the angles between spheroid directions and {\mbox{\boldmathr}}_{ij}={\mbox{\boldmathr}}_{i}-{\mbox{\boldmathr}}_{j} as EPL ; Takae-ori

[TABLE]

where $i\in 1$ in the first term, $j\in 1$ in the second term, and $\eta$ represents the molecular anisotropy. For $\eta>0$, we have $0\leq A_{ij}\leq 2\eta$, which vanishes for {\mbox{\boldmathn}}_{i}\cdot{\mbox{\boldmathr}}_{ij}=0 for $j\in 2$ (and {\mbox{\boldmathn}}_{j}\cdot{\mbox{\boldmathr}}_{ij}=0 for $j\in 1$). We assume a relatively small size difference and mild anisotropy as

[TABLE]

Then, at density $0.84\sigma_{1}^{-3}$, our system forms a crystal without phase separation and isotropic-nematic phase transition Latz ; EPL . For larger $\sigma_{2}/\sigma_{1}$ and $\eta$, the latter processes may take place during slow quenching from liquid. Because $U_{\rm LJ}$ is minimized at $r_{ij}=(2+2A_{ij})^{1/6}\sigma_{\alpha\beta}$ for fixed {\mbox{\boldmathn}}_{i} and {\mbox{\boldmathn}}_{j}, we regard the anisotropic particles as spheroids with aspect ratio $(1+2\eta)^{1/6}=1.24$. Notice that our potential is similar to the Gay-Berne potential for rodlike molecules Berne .

The second term $U_{\rm w}$ in Eq.(2) is the sum of strongly repulsive, wall potentials as apply

[TABLE]

We set $w=e^{40}\epsilon$ and $\xi_{\rm w}=0.01\sigma_{1}$ to make the potentials hardcore-like. Then, the distances between the dipole centers and the walls become longer than $0.4\sigma_{1}$.

II.2 Electrostatic energy and canonical distribution

We assume permanent dipolar moments {\mbox{\boldmath\mu}}_{i} along the spheroid direction {\mbox{\boldmathn}}_{i} ($1\leq i\leq N_{1}$) written as

[TABLE]

where $\mu_{0}$ is a constant dipole moment. There is no induced dipole moment. The electric potential \Phi({\mbox{\boldmathr}}) can be defined away from the dipole positions {\mbox{\boldmathr}}\neq{\mbox{\boldmathr}}_{i}. We impose the metallic boundary condition at $z=0$ and $H$:

[TABLE]

where $\Delta\Phi$ is the applied potential difference and $E_{\rm a}=\Delta\Phi/H$ is the applied electric field. In this paper, we perform simulation by controlling $\Delta\Phi$ (or $E_{\rm a}$). In our scheme, $\Delta\Phi$ can be nonstationary.

The boundary condition (8) is realized by the surface charge densities at $z=0$ and $H$ (see Appendix A). As a mathematical convenience, we instead introduce image dipoles outside the cell for each dipole {\mbox{\boldmath\mu}}_{i} at {\mbox{\boldmathr}}_{i}=(x_{j},y_{j},z_{j}) in the cell. As in Fig.1(a), we first consider those at {\mbox{\boldmathr}}_{i}-2Hm_{z}{\mbox{\boldmathe}}_{z} ($m_{z}=\pm 1,\pm 2,\cdots$) with the same moment {\mbox{\boldmath\mu}}_{i}, where {\mbox{\boldmathe}}_{z} is the unit vector along the $z$ axis. Second, at {\bar{\mbox{\boldmathr}}}_{i}-2Hm_{z}{\mbox{\boldmathe}}_{z} $(m_{z}=0,\pm 1,\pm 2,\cdots)$, we consider those with the image moment given by

[TABLE]

where {\bar{\mbox{\boldmathr}}}_{i}=(x_{i},y_{i},-z_{i}) is the image position closest to the bottom wall. For {\mbox{\boldmathr}}\neq{\mbox{\boldmathr}}_{i}, the real and image dipoles and the applied field yield the following potential,

[TABLE]

where {\mbox{\boldmathr}}\neq{\mbox{\boldmathr}}_{i}, {\mbox{\boldmathg}}({\mbox{\boldmathr}})=r^{-3}{\mbox{\boldmathr}}, and {\mbox{\boldmathh}}=(Lm_{x},Lm_{y},2Hm_{z}) with $m_{x},m_{y}$, and $m_{z}$ being integers. Here, the first term is periodic in three dimensions (3D). Along the $z$ axis the period is $2H$ because of the summation over $m_{z}$ or over the image dipoles. We confirm that the first term in Eq.(10) vanishes at $z=0$ and $H$ with the aid of Eq.(9).

At fixed $E_{\rm a}$, the total electrostatic energy $U_{\rm d}$ in Eq.(2) is now written in terms of {\mbox{\boldmathr}}_{i} and {\mbox{\boldmath\mu}}_{i} asKlapp ; apply

[TABLE]

Here, {\stackrel{{\scriptstyle\leftrightarrow}}{{\mbox{\boldmath\cal T}}}}({\mbox{\boldmathr}}) is the dipolar tensor with its $\alpha\beta$ component being $\delta_{\alpha\beta}/r^{3}-3x_{\alpha}x_{\beta}/r^{5}$. In the first term, the self-interaction contributions ({\mbox{\boldmathh}}={\mbox{\boldmath0}} and $i=j$) are removed in $\sum_{i\in 1,j\in 1}^{\prime}$. In the second term, we set {\bar{\mbox{\boldmathr}}}_{ij}={\mbox{\boldmathr}}_{i}-{\bar{\mbox{\boldmathr}}}_{j}=(x_{i}-x_{j},y_{i}-y_{j},z_{i}+z_{j}). In the last term, $M_{z}$ is the $z$ component of the total polarization,

[TABLE]

For each dipole $i$, the electrostatic force is given by {\mbox{\boldmathF}}_{i}^{\rm e}=-\partial U/\partial{\mbox{\boldmathr}}_{i} and the local electric field by

[TABLE]

We can also obtain {\mbox{\boldmathE}}_{i} by subtracting the self contribution {\mbox{\boldmathg}}({\mbox{\boldmathr}}-{\mbox{\boldmathr}}_{i})\cdot{\mbox{\boldmath\mu}}_{i} from \Phi({\mbox{\boldmathr}}) in Eq.(10) as

[TABLE]

We consider the Hamiltonian ${\cal H}={\cal K}+U$, where ${\cal K}$ is the total kinetic energy. In our model, the applied field $E_{\rm a}$ appears linearly in $U_{\rm d}$ in Eq.(11). Then, we find

[TABLE]

where ${\cal H}_{0}$ is the Hamiltonian for $E_{\rm a}=0$. This form was assumed in the original linear response theoryKubo . For stationary $E_{\rm a}$, the equilibrium average, denoted by ${\langle{\cdots}\rangle}_{\rm e}$, is over the canonical distribution $\propto\exp(-{\cal H}/k_{\rm B}T)$. Then, for any variable ${\cal A}$ (independent of $E_{\rm a}$), its equilibrium average $\langle{{\cal A}}\rangle_{\rm e}$ changes as a function of $E_{\rm a}$ asOnukibook

[TABLE]

where $T$ is fixed in the derivative and $\delta M_{z}=M_{z}-\langle{M_{z}}\rangle_{\rm e}$. For the average polarization $P_{z}=\langle{M_{z}}\rangle_{\rm e}/V$, we consider the differential susceptibility $\chi_{\rm dif}=\partial P_{z}/\partial E_{\rm a}$. In equilibrium, it is related to the variance of $\delta M_{z}$ as

[TABLE]

As $E_{\rm a}\to 0$, $\chi_{\rm dif}$ tends to the susceptibility $\chi=(\varepsilon-1)/4\pi$ 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 $T$, but we do not obtain Eq.(17) at low $T$ 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 $\cal K$ depends on the translational velocities {\dot{\mbox{\boldmathr}}}_{i}={d{\mbox{\boldmathr}}_{i}}/{dt} ($i=1,\cdots,N$) and the angular velocities {\dot{\mbox{\boldmathn}}}_{i}={d{\mbox{\boldmathn}}_{i}}/{dt} ($i=1,\cdots,N_{1}$) as

[TABLE]

where $m$ is the mass common to the two species, and $I_{1}$ is the moment of inertia. We set $I_{1}=0.125m\sigma_{1}^{2}$ in this paper. The Newton equations for {\mbox{\boldmathr}}_{i} are given by

[TABLE]

where {\ddot{\mbox{\boldmathr}}}_{i}={d^{2}{\mbox{\boldmathr}}_{i}}/{dt^{2}}. On the other hand, the Newton equations for {\mbox{\boldmathn}}_{i} ($1\leq i\leq N_{1}$) are of the formapply ; EPL ,

[TABLE]

where {\ddot{\mbox{\boldmathn}}}_{i}={d^{2}{\mbox{\boldmathn}}_{i}}/{dt^{2}}, \stackrel{{\scriptstyle\leftrightarrow}}{{\mbox{\boldmathI}}} is the unit tensor, and {\mbox{\boldmathE}}_{i}^{\rm eff}=-\partial U/\partial{\mbox{\boldmath\mu}}_{i} is the local orientating field on dipole $i$. The left hand side of Eq.(20) is perpendicular to {\mbox{\boldmathn}}_{i} from {\mbox{\boldmathn}}_{i}\cdot{\ddot{\mbox{\boldmathn}}}_{i}+|{\dot{\mbox{\boldmathn}}}_{i}|^{2}=0. The right hand side vanishes if {\mbox{\boldmathE}}_{i}^{\rm eff} is parallel to {\mbox{\boldmathn}}_{i}. From Eqs.(19) and (20) the Hamiltonian ${\cal H}={\cal K}+U$ changes as ${d}{\cal H}/dt=-M_{z}{d}E_{\rm a}/dt$ (without thermostats). Thus, $\cal H$ is conserved for stationary $E_{\rm a}$.

At low $T$, we have \mu_{0}|{\mbox{\boldmathE}}_{i}^{\rm eff}|\gg k_{\rm B}T for most $i\in 1$, where {\mbox{\boldmathn}}_{i} is nearly parallel to {\mbox{\boldmathE}}_{i}^{\rm eff}. From Eq.(2) we set

[TABLE]

where {\mbox{\boldmathE}}_{i} is the long-range dipolar part in Eq.(13) and {\mbox{\boldmathE}}_{i}^{\rm ste} is the short-range steric part from the orientation dependence of $U_{\rm LJ}$ in Eq.(3). Some calculations give

[TABLE]

where main contributions arise from neighbors $j$ with $r_{ij}\lesssim\sigma_{\alpha\beta}$. These neighbor impurities ($j\in 2$) yield local random pinning fields (see Fig.3(a)).

II.4 Simulation method

We integrated Eqs.(19) and (20) for $N=N_{1}+N_{2}=8000$. We used the 3D Ewald method on the basis of $U_{\rm d}$ in Eq.(11) Hautman ; Klapp ; apply ; Takae1 ; Takae2 . To realize crystal, we slowly cooled the system from a liquid above the melting temperature ($\sim\epsilon/k_{\rm B}$) at density $N/V=0.84\sigma_{1}^{-3}$. In crystal, there is no translational diffusion. We attached Nosé-Hoover thermostats to the particles in the layer regions $z<2\sigma_{1}$ and $H-z<2\sigma_{1}$. We fixed the cell volume at $V=L^{2}H$ with $H=L=21.2\sigma_{1}$ mostly, but we slightly varied $H$ in time to obtain the field-induced strain in Sec.V.

In our system, the dielectric response strongly depends on the dipole moment $\mu_{0}$ in Eq.(7)Johnson1 , so we present our results for $\mu_{0}=0.8$ and 1.6 in units of $(\epsilon\sigma_{1}^{3})^{1/2}$. For exampleJohnson0 ; Johnson1 , if $\epsilon/k_{\rm B}=100$ K and $\sigma_{1}=5~{}{\rm\AA}$, these values of $\mu_{0}$ are $1.05$ D and 2.10 D, respectively.

We measure space and time in units of $\sigma_{1}$ and

[TABLE]

Units of $T$, electric potential, and electric field are $\epsilon/k_{\rm B}$, $(\epsilon/\sigma_{1})^{1/2}$, and $(\epsilon/\sigma_{1}^{3})^{1/2}$, respectively. For $\epsilon/k_{\rm B}=100$ K and $\sigma_{1}=5~{}{\rm\AA}$, we have $(\epsilon/\sigma_{1})^{1/2}=0.16$ V, $(\epsilon/\sigma_{1}^{3})^{1/2}=0.32$ V$/$nm, and $e=18.3(\epsilon\sigma_{1})^{1/2}$ (elementary charge).

Because of heavy calculations of electrostatics we performed a single simulation run for each parameter set. Then, $\langle{\cdots}\rangle$ denotes the time average (not the ensemble one). We also do not treat slow aging processesaging ; aging0 ; Kob , for which very long simulation time is needed.

III Ferroelectric transition for $c=0$

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

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

[TABLE]

where $P_{1}(x)=x$ and $P_{2}(x)=(3x^{2}-1)/2$. We sum over neighbor dipoles $j$ with $r_{ij}<1.4$, where $Z_{i}$ is their number. Then, $Q_{1i}$ represents the local dipolar order and $Q_{2i}$ the local quadrupolar orderBinder ; ori . These variables will be used also for ferroelectric glass with $c>0$.

In Fig.2, we examine the transition by slowly lowering $T$ for ${c=0}$, $\mu_{0}=0.8$, and $\Delta\Phi=0$. In (a), we plot the averages $\langle{Q_{1}}\rangle=\sum_{i\in 1}Q_{1i}/N_{1}$ and $\langle{Q_{2}}\rangle=\sum_{i\in 1}Q_{2i}/N_{1}$. Here, the transition is steep but gradual due to the finite-size effect imposed by the metal walls at $z=0$ and $H$. In our system, the spheroidal particles form a fcc crystal in the plastic crystal phase Binder ; ori ; EPL in the range $0.3<T<1$. For lower $T$, a polycrystal with eight rhombohedral variants appears, where the spheroid directions {\mbox{\boldmathn}}_{i} are along $\langle{111}\rangle$ except those near the interfaces.

In the transition range $0.22\lesssim T\lesssim 0.30$, the system is composed of disordered and ordered regions with sharp interfaces. We give snapshots of relatively ordered regions with $Q_{1i}>0.6$ at (b) $T=0.28$ and (c) $T=0.25$, where we pick up (b) 10$\%$ and (c) $30\%$ of the total dipoles. These patterns are stationary in our simulation time intervals. In (d), at $T=0.05$, we give a snapshot of polycrystal state with eight variants.

The rhombohedral structure is characterized by the angles $\pi/2\pm\alpha$ of its lozenge faces of a unit cell. At low $T$, we find $\alpha\cong 5^{\circ}$ for $\mu_{0}=0.8$ but $\alpha\cong 1^{\circ}$ for $\mu_{0}=1.6$. See Sec.VA for the reason of this $\mu_{0}$ dependence.

IV Ferroelectric transition for $c>0$

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 Takae-ori ; EPL , 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 Jun .

In Fig.3(a), we display the dipole directions for c=0.2 on a (111) plane at $z=0.75H$. Many of them tend to align in the directions parallel to the impurity surfaces or perpendicular to {\mbox{\boldmathr}}_{ij} ($j\in 2$), because $A_{ij}=0$ for {\mbox{\boldmathn}}_{i}\cdot{\mbox{\boldmathr}}_{ij}=0 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 Ann ; Bokov2 ; Bokov4 .

In Fig.3(b), we display all the impurities in the cell for $c=0.2$, 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 $(\gtrsim 10)$ 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 effortWel ; Wel1 ; Ban has been made to determine the distribution of the B-site ions (Mg2+ and Nb5+ for PMN) using effective atom-atom interactions, while Burton et al. Burton1 ; Burton2 ; Burton3 demonstrated strong influence of compositional heterogeneity on the PNRs.

IV.2 Diffuse transition toward ferroelectric glass

The dipole moment $\mu_{0}$ determines relative importance of the dipolar and steric parts, {\mbox{\boldmathE}}_{i} and {\mbox{\boldmathE}}_{i}^{\rm ste}, in the orientating field in Eq.(21), since they are proportional to $\mu_{0}$ and $\mu_{0}^{-1}$, respectively, for $E_{\rm a}=0$. For example, for $c=0.2$ and $\Delta\Phi=0$, we average ({|{\mbox{\boldmathE}}_{i}|}, {|{\mbox{\boldmathE}}_{i}^{\rm ste}|}) over all the dipoles to obtain ($2.0$, $5.7$) for $\mu_{0}=0.8$ and $T=0.1$ and ($4.6$, $3.4$) for $\mu_{0}=1.6$ and $T=0.2$. Thus, {\mbox{\boldmathE}}_{i} is more important for larger $\mu_{0}$ in the dipole orientations. Here, the amplitude of the local electric field |{\mbox{\boldmathE}}_{i}| is mostly of order $4\pi\mu_{0}n_{1}/3=2.8\mu_{0}$, where $n_{1}=(1-c)N/V$. This large size of {\mbox{\boldmathE}}_{i} is realized within mesoscopic PNRs.

In Fig.4, we examine the transition with $c=0.2$ and $\Delta\Phi=0$ for the two cases $\mu_{0}=0.8$ and $1.6$, where the net polarization nearly vanishes. At each $T$, we waited for a time $\Delta t\sim 5\times 10^{4}$. In (a) and (b) we show gradual $T$ dependence of $\langle{Q_{\ell}}\rangle=\sum_{i\in 1}Q_{\ell i}(t_{0})/N_{1}$. They take appreciable values in the presence of small PNRs. In (c) and (d), we also show their variances,

[TABLE]

The orientation fluctuations are frozen at large sizes at low $T$. We also see that $\langle{Q_{1}}\rangle$ in (a) and $\langle{(\delta Q_{1})^{2}}\rangle$ in (c) exhibit small maxima at low $T$, but they should disappear in the ensemble averages.

In Fig.5, we plot the time-correlation functions $C_{1}(t)$ for one-body angle changes defined by

[TABLE]

where the average is taken over the initial time $t_{0}$. In (a) and (b), the angle changes slow down with lowering $T$. We define the reorientation time $\tau_{1}$ by

[TABLE]

where 0.1 is smaller than the usual choice $e^{-1}$ since $C_{1}(t)$ decays considerably in the initial thermal stage for not very low $T$. The PNRs are broken on this timescale. In (c) and (d), we display $\tau_{1}$ vs $1/T$. where $\tau_{1}$ can well be fitted to the Vogel-Fulcher form Kob ,

[TABLE]

Here, $\tau_{10}$, $T_{1}$, and $D_{1}$ are constants with $(D_{1},T_{1})$ being $(0.89,0.19)$ for $\mu_{0}=0.8$ and $(4.7,0.32)$ for $\mu_{0}=1.6$.

IV.3 Dielectric permittivity

We next examine the dielectric permittivity. We calculated its real part $\varepsilon^{\prime}(\omega,T)$ and imaginary part $\varepsilon^{\prime\prime}(\omega,T)$ as functions of $T$ and the frequency $\omega$ by applying small ac field in the linear response regime (see Appendix B).

In Fig.6, we show $\varepsilon^{\prime}$ and the ratio $\varepsilon^{\prime\prime}/\varepsilon^{\prime}$ vs $T$ at several low frequencies for the two cases $\mu_{0}=0.8$ (left) and 1.6 (right). In (a) and (b), $\varepsilon^{\prime}$ increases with decreasing $\omega$ and exhibit a broad maximum at a temperature $T=T_{\rm m}(\omega)$ for each $\omega$. With decreasing $\omega$, $T_{\rm m}(\omega)$ decreases (with weaker dependence for smaller $\omega$) and the peak height $\varepsilon^{\prime}_{\rm m}(\omega)=\varepsilon^{\prime}(\omega,T_{\rm m}(\omega))$ increases. For $T>T_{\rm m}$, we have $\omega<\tau_{1}^{-1}$, so $\varepsilon^{\prime}(\omega,T)$ tends to the linear dielectric constant $\varepsilon(T)$. However, for $T<T_{\rm m}$, $\varepsilon^{\prime}$ decreases to zero with lowering $T$ or increasing $\omega$, where the response of the PNRs to small ac field decreases. On the other hand, $\varepsilon^{\prime\prime}/\varepsilon^{\prime}$ exhibits a maximum for each $\omega$ and shifts to a lower temperature with lowering $\omega$. These behaviors characterize ferroelectric glass Cr ; Cr1 ; Smo ; Klee ; Klee1 ; Zhang ; Cowley ; Bokov ; Bokov3 ; Samara ; Loidl3 ; Ho2 . Similar behaviors were found for the frequency-dependent magnetic susceptibilities in spin glassMydosh . Furthermore, in Appendix B, we will present analysis of $\varepsilon^{\prime}$ and $\varepsilon^{\prime\prime}$ for $\omega\gtrsim\tau_{1}^{-1}$ at relatively high $T$ on the basis of the linear response theoryKubo .

We write the inverse relation of $T=T_{\rm m}(\omega)$ as

[TABLE]

leading to $\varepsilon^{\prime}(\omega_{\rm m},T)=\varepsilon_{\rm m}^{\prime}(\omega_{\rm m})$. Here, $T/T_{\rm m}>1$ ($<1$) holds for $\omega/\omega_{\rm m}<1$ ($>1$). In (c) and (d) of Fig.5, we compare the inverse $1/\omega_{\rm m}(T)$ and $\tau_{1}(T)$ in Eq.(27) for $\mu_{0}=0.8$ and 1.6. We find $\omega_{\rm m}\sim\tau_{1}^{-1}$. Thus, $\omega_{\rm m}(T)$ represents a characteristic frequency of the dipole reorientations. Previously, for relaxors and spin glasses, Stringer et al.Ra nicely fitted $1/\omega_{\rm m}(T)$ to the Vogel-Fulcher form, which is in accord with (c) and (d) of Fig.5.

Our system is ergodic at relatively high $T$, but becomes nonergodic as $T$ is lowered. The boundary between these two regimes weakly depends on the observation time. In Fig.7(a), $M_{z}(t)$ evolves on a wide range of time scales in a time interval with width $6\times 10^{4}$ for $c=0.2$, $\Delta\Phi=0$, and $\mu_{0}=1.6$. At $T=0.5$, its time average becomes small, but its fluctuations are large. In contrast, at $T=0.35$, it remains negative around $-800=-14(Vk_{\rm B}T)^{1/2}$, on which smaller thermal fluctuations with faster time scales are superimposed. Note that the ensemble average of $M_{z}$ should vanish at any $T$ for $\Delta\Phi=0$.

For a single simulation run, we consider the time average of the normalized polarization variance, written as $\chi_{\rm fl}$. To avoid confusion, we define it explicitly as

[TABLE]

We set ${\langle{{\cal A}}\rangle}_{\rm time}=\int_{t_{1}}^{t_{2}}dt{\cal A}(t)/\Delta t$ with $\Delta t=t_{2}-t_{1}(\sim 5\times 10^{4}$ here) for any time-dependent variable ${\cal A}(t)$. This averaging procedure has already been taken for the quantities in Figs.4 and 5. In the nonergodic $T$ range, ${\langle{M_{z}}\rangle}_{\rm time}$ remains nonvanishing even for $\Delta\Phi=0$, while $\chi_{\rm fl}$ arises from the (thermal) dynamical fluctuations and tends to zero as $T\to 0$. In Fig.7(b), we plot numerical results of $\varepsilon_{\rm fl}=1+4\pi\chi_{\rm fl}$ for $\Delta\Phi=0$ and $\varepsilon^{\prime}=1+4\pi\chi^{\prime}$ at $\omega=2\pi\times 10^{-4}$ as functions of $T$. These two curves nearly coincide for $T>T_{\rm m}$ yielding $\varepsilon(T)$, but $\varepsilon_{\rm fl}$ is considerably larger than $\varepsilon^{\prime}$ for $T<T_{\rm m}$. In their simulation, Burton et al.Burton1 ; Burton2 calculated a dielectric constant from polarization fluctuations, which corresponds to $\varepsilon_{\rm fl}$ in our case.

In Fig.7(b), $\varepsilon_{\rm fl}$ and $\varepsilon^{\prime}$ steeply grow as $T\to T_{\rm m}$. From the curves of $1/\chi^{\prime}$ and $1/\chi_{\rm fl}$ in its inset, $\chi^{\prime}=(\varepsilon^{\prime}-1)/4\pi$ and $\chi_{\rm fl}$ can fairly be fitted to the Curie-Weiss form,

[TABLE]

with $A_{0}\cong 1.2$ and $T_{0}\cong 0.47(\cong T_{\rm m}$ at $\omega=2\pi\times 10^{-4}$) for $T\gtrsim 0.55$. At $T=0.5$, however, we find $\chi^{\prime}\cong 13$ and $\chi_{\rm fl}\cong 23$. In experiments, the behavior (31) was found for orientational glassori , but a marked deviation was detected close to $T_{\rm m}$ for relaxorsSamara ; Bokov ; Cr1 . Thus, if $T$ is somewhat above $T_{0}$, our polarization fluctuations resemble the critical fluctuations in systems near their critical point Cowley ; Samara . In our disordered system, these near-critical fluctuations are slowed down and eventually frozen as $T$ is further lowered, as in relaxors. This can also be seen in (c) and (d) of Fig.4. Furthermore, for $T\lesssim T_{0}$, there is a tendency of interface formation between adjacent PNRs for $c\lesssim 0.2$, which will be discussed in future.

For relaxors, Stock et al.Stock divided the scattering intensity into frozen and dynamic parts, where the former (latter) increases (decreases) with lowering $T$. Similar arguments of nonergodicity were made for polymer gelsPusey ; Matsuo , 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 Matsuo ; Onukibook .

We next confirm Eq.(17) by increasing $E_{\rm a}$ at $T=0.6$ with $c=0.2$ and $\mu_{0}=1.6$, where the observation time is much longer than $\tau_{1}\sim 60$. In Fig.7(c), we compare the differential formula $\varepsilon_{\rm dif}=1+4\pi dP_{z}/dE_{\rm a}$ and the fluctuation formula $\varepsilon_{\rm fl}=1+4\pi\chi_{\rm fl}$ 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 $T$, the polarization fluctuations are suppressed with increasing $E_{\rm a}$.

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 $T$, they have finite lifetimes (within observation times)Bokov ; Bl . This feature is illustrated in two snapshots in Fig.8, which were taken at two times separated by $10^{4}$ in the same simulation run. They display the dipoles with $Q_{1i}>0.6$ for $T=0.22$, $c=0.2$, and $\mu_{0}=0.8$. These two patterns are very different, so their lifetime is shorter than $10^{4}$. In fact, $\tau_{1}$ is of order $10^{3}$ at $T=0.22$ in Fig.5(c).

The PNRs are frozen with lowering $T$. In Fig.9, we give examples for $c=0.1$, 0.2, and 0.3 with $T=0.05$, $\mu_{0}=0.8$, and $\Delta\Phi=0$. The left panels display the particles on the boundaries ($x,y,{\rm or}~{}z=L$), while the right ones the relatively ordered dipoles with $Q_{1i}>0.7$. The dipoles depicted in the latter amount to 37, 20, and 13$\%$ of the total dipoles from above. For $c=0.1$, we can see well-defined ordered domains consisting of eight variants, whose interfaces are trapped at impuritiesEPL ; Takae-ori (see Fig.10(a)). These domains are broken up into smaller PNRs with increasing $c$. For $c=0.2$, the PNRs mostly take compressed, plate-like shapes under the constraint of the spatially correlated impurities (see Fig.10 also). For $c=0.3$, 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 $i$ satisfies $Q_{1i}>0.7$ and $r_{ij}<1.4$ for some $j$ within the same PNR. In Fig.9, the dipole number in a PNR is $2200,270$, and $56$ on the average from above. Thus, the connectivity of the PNRs sensitively depends on $c$. In the following, we treat the case $c=0.2$.

IV.5 Single polar nanoregion and local electric field

We visualize individual PNRs frozen at low $T$. When the system is composed of PNRs, the local electric field {\mbox{\boldmathE}}_{i} in Eq.(13) arises mainly from the dipoles within the same PNR in the bulk. Its amplitude is of order $4\pi\mu_{0}n_{1}/3\cong 2.8\mu_{0}$ for not very large $\Delta\Phi$.

In Fig.10, we pick up (a) a single PNR for $\Delta\Phi=0$ at the cell center and (b) another one for $\Delta\Phi=3$ ($E_{\rm a}=0.14$) near the upper wall, where $c=0.2$ and $T=0.05$. 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) $(120,70)$ and (b) $(70,50)$ 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 {\mbox{\boldmathE}}_{i} and {\mbox{\boldmath\mu}}_{i}, where they are nearly parallel. Here, we divide the dipolar part of {\mbox{\boldmathE}}_{i} into the contributions from those inside and outside the PNR, written as {\mbox{\boldmathE}}_{i}^{\rm d}(\rm{in}) and {\mbox{\boldmathE}}_{i}^{\rm d}(\rm{out}). Then, Eq.(A9) in Appendix A gives

[TABLE]

where the last term arises from the surface charges. In (a), we find {\mbox{\boldmathE}}_{i}\cong{\mbox{\boldmathE}}_{i}^{\rm d}(\rm{in}), which occurs mostly for the dipoles in the interior of PNRs in the bulk. In (b), on the other hand, we find {\mbox{\boldmathE}}_{i}\cong{\mbox{\boldmathE}}_{i}^{\rm sur}, where |{\mbox{\boldmathE}}_{i}^{\rm sur}| is of the same order as $4\pi{\bar{\sigma}}_{0}=2.4$ and is much larger than $E_{\rm a}=0.14$. Here, ${\bar{\sigma}}_{0}$ is the mean surface charge density at $z=0$. For example, if we set $\epsilon/k_{\rm B}=100$ K and $\sigma_{1}=5~{}{\rm\AA}$, we have ${\bar{\sigma}}_{0}\sim 0.01e/\sigma_{1}^{2}$. for $\Delta\Phi=3$.

IV.6 Orientation near metal surface

As can be seen in Figs.9 and 10(b), the dipoles next to the walls are parallel or antiparallel to the $z$ axis (along $[111]$), whose distances from the walls are about 0.5. This is due to their interaction with the image dipoles in the walls (see Appendix A)Takae1 ; Takae2 ; apply . For $\Delta\Phi=0$, these two orientations appear equally on the average due to the top-tail symmetry of our spheroidal dipoles. For $\Delta\Phi\neq 0$, one of them is more preferred than the other. In Fig.11, we show the particles in the first and second $(111)$ layers in applied field with $\Delta\Phi=3$, where $c=0.2$, $\mu_{0}=0.8$, and $T=0.05$. The parallel and antiparallel orientations appear in the first layer, but the other oblique ones also appear in the second layer. We shall see that the corresponding surface charge density $\sigma_{0}(x,y)$ at $z=0$ is highly heterogeneous in Fig.17 in Appendix A.

In our crystal case, the $\ell$-th layer is given by $\ell-1<z<\ell$, since the separation between adjacent $(111)$ planes is close to 1. Here, we consider the average of $n_{zi}=\cos\theta_{i}$ over the dipoles in the $\ell$-th layer and write it as $\langle{n_{z}}\rangle_{\ell}$. In Fig.11, it is 0.30 for $\ell=1$, $0.34$ for $\ell=2$, and $\langle{n_{z}}\rangle_{\rm b}=0.33$ for $\ell\gg 1$. These values are close, so the surface effect on the polarization is weak in this case of our model. The excess potential drop near the bottom wall is given by $4\pi\mu_{0}n_{1}\sigma_{1}\sum_{\ell}[\langle{n_{z}}\rangle_{\rm b}-\langle{n_{z}}\rangle_{\ell}]\cong 0.13$, which is much smaller than the total drop $\Delta\Phi=3$. In contrast, for highly polar liquids such as water Hautman ; Takae2 ; Willard , a significant potential drop appears in a microscopic (Stern) layer on a solid surface even without ion adsorption.

V Polarization and strain in applied electric field

V.1 Applying electric field along $[111]$ at fixed stress

In this section, we give results of cyclic changes of $\Delta\Phi=HE_{\rm a}$ for $c=0.2$ and $\mu_{0}=1.6$. We also calculated the response with $\mu_{0}=0.8$ (not shown here). For these two $\mu_{0}$ values, the characteristic features are nearly the same, but the response sizes are very different. That is, the dielectric response for $\mu_{0}=1.6$ is larger than that for $\mu_{0}=0.8$ by one order of magnitude as in Fig.6, while the field-induced strain for $\mu_{0}=1.6$ is about $20\%$ of that for $\mu_{0}=0.8$. See the last paragraph of Sec.III for the rhombohedral angles in our simulation. Using a barostat, we fixed the $zz$ component of the average stress and varied the cell width $H(t)$ to calculate the field-induced strain. The lateral cell length was fixed at $L$.

In our model, dipole alignment along $[111]$ yields both steric repulsion and dipolar attraction between adjacent $(111)$ planes. Their relative importance depends on $\mu_{0}$. If the former is larger (smaller) than the latter, an expansion (a shrinkage) of the cell width $H$ occurs for $\Delta\Phi\neq 0$. Note that the dipolar interaction between two dipoles at {\mbox{\boldmathr}}_{i} and {\mbox{\boldmathr}}_{j} aligned along the $z$ axis is attractive (repulsive) if the angle between their relative vector {\mbox{\boldmathr}}_{i}-{\mbox{\boldmathr}}_{j} and the $z$ axis is smaller (larger) than $\cos^{-1}(1/\sqrt{3})$.

We increased $\Delta\Phi(t)$ from 0 to 10, decreased to $-10$, and then increased again to 10 at fixed $T$ without dislocation formation. The changing rate ${\dot{\Phi}}=d(\Delta\Phi)/dt$ was $\pm 1.5\times 10^{-3}$ $\epsilon/(\sigma_{1}^{3}m)^{1/2}$. The average pressure along the $z$ axis was 0.4 at $T=0.2$ and 3.6 at $T=0.6$ in units of $\epsilon/\sigma_{1}^{3}$, while the lateral one increased by 0.8 for a change of $\Delta\Phi$ from 0 to 10. The $H(t)$ changed from $H(0)=L$ at most by $2\%$. We calculated the average polarization and strain for $t>0$ given by

[TABLE]

We also calculated the mean surface charge density ${\bar{\sigma}}_{0}$ at $z=0$ to confirm Eq.(A5) in Appendix A (see Fig.12).

V.2 Polarization response

In Fig.12, we plot $4\pi{P}_{z}$ vs $\Delta\Phi$ for (a) $T=0.6$ and (b) $T=0.2$. At $\Delta\Phi=10$, we have $P_{z}\sim 10$ and $\langle{n_{z}}\rangle=P_{z}/n_{1}\mu_{0}\sim 0.75$ in (a) and (b), where $n_{zi}$ is the $z$ component of {\mbox{\boldmathn}}_{i} and $n_{1}=N_{1}/V$ is $0.672$. In (a), there is no hysteresis and the initial slope yields $\varepsilon_{\rm dif}=1+4\pi dP_{z}/dE_{\rm a}\cong 153$. In (b), marked hysteresis appears, where $\varepsilon_{\rm dif}\cong 13$ initially at point I, but is about $150$ between two points A and B ($\Delta\Phi:4\to 6$). Here, the initial $P_{z}$ at I is slightly negative as a frozen fluctuation (see the curve at $T=0.35$ in Fig.7(a)). At point C on the vertical axis we have a remnant polarization $P_{\rm R}=0.72$ with $\Delta\Phi=0$. For any $T$, $\varepsilon_{\rm dif}$ from the initial slope at $E_{\rm a}=0$ nearly coincides with $\varepsilon^{\prime}$ at $\omega=2\pi\times 10^{-4}$ in Fig.7(b) (equal to 50 at $T=0.4$ and to 10 at $T=0.1$). The curves in (a) and (b) closely resemble those in various ferroelectric systems Bl ; Mori ; Furukawa ; Zhang ; piezo .

The field-induced change from A to B in (b) is very steep with large $\chi_{\rm dif}$. In (c), we thus display the dipoles with large angle changes: $\Delta n_{zi}=n_{zi}(t_{B})-n_{zi}(t_{A})>1$, where $t$ is $t_{A}$ at A and $t_{B}$ at B. Collective reorientations are marked in this time interval. In (d), for three intervals, we plot the distribution function $P_{\rm a}(u)=\sum_{i\in 1}{\delta(u-\Delta n_{zi})}/N_{1}$ for $\Delta n_{zi}=n_{zi}(t_{1})-n_{zi}(t_{0})$, where we use an appropriately smoothed $\delta$-function. Small angle changes are dominant in the first interval $t_{I}<t<t_{A}$ (where $t=t_{I}$ at I), but large angle changes are dominant in the next interval $t_{A}<t<t_{B}$.

In (b), the initial point I (at $t=t_{I}$) of the cycle represents an arrested state with frozen fluctuations realized by zero-field cooling. It can no longer be reached once a large field is applied. The corresponding states have been realized in many systems (see Sec.VI). In the two states at I and C, the polarization directions are very different, but the values of the potential energy $U$ in Eq.(2) are close as $-8.07N\epsilon$ at I and $-8.04N\epsilon$ at C. We can also see that the quadrupolar order parameters $Q_{2i}(t)$ in Eq.(24) do not change much for most $i$ during the cycle despite large changes in {\mbox{\boldmathn}}_{i}. For example, the mean square difference $\sum_{i}[Q_{2i}(t_{I})-Q_{2i}(t)]^{2}/N_{1}$ for time interval $[t_{I},t]$ is $0.047$, $0.103$, and $0.086$ at $t=t_{A}$, $t_{B}$, and $t_{C}$ (which are the times at A, B, and C), respectively, where the variance $\sum_{i\in 1}[\delta Q_{2i}(t)]^{2}/N_{1}$ for $\delta Q_{2i}=Q_{2i}-\langle{Q_{2}}\rangle$ remains of order 0.08 (see Eq.(25) and Fig.4d)

Between A and B in (b), we found an increase in the polarization variance $\langle{(\delta M_{z})^{2}}\rangle(t)$ $(t_{A}\lesssim t\lesssim t_{B}$). For relaxors, Xu et al. Shi2 detected an increase in the diffuse scattering in the field range with large $\chi_{\rm dif}$. We should then examine the scattering amplitude between A and B. In addition, when $\Delta\Phi$ was held fixed at 4.0 (at A), we observed slow reorientations leading to coarsening of PNRsaging ; aging0 . These effects will be studied in future.

V.3 Field-induced strain

In our model, the heterogeneity in the strain is marked because of dilation of PNRs along $\langle{111}\rangle$, though it is milder than that of the polarization. To illustrate this effect, we define a local strain $\gamma_{i}$ along the $z$ axis for each particle $i$ (including the impurities) by

[TABLE]

where the summation is over other $j$ with $r_{ij}<1.4$ and $|z_{j}-z_{i}|>0.2$, ${W}_{i}$ is the number of these neighbors, and $a_{111}(\cong 1.0)$ is the average spacing between two consecutive $(111)$ planes. From these conditions, the $(111)$ plane containing $j$ is adjacent to that containing $i$. The particle average $\sum_{i}\gamma_{i}/N$ nearly coincides with $\bar{\gamma}$ in Eq.(33).

In Fig.13, we plot $\bar{\gamma}$ vs $\Delta\Phi$ with $\mu_{0}=1.6$ in the same simulation run as in Fig.12. We find (a) a cusp curve at $T=0.6$ and (b) a butterfly-like curve at $T=0.2$. In (b), $\bar{\gamma}$ becomes slightly negative at $\Delta\Phi\cong\pm 4$. These two curves resemble those in the previous experiments Furu1 ; Zhang ; piezo . In (c), we pick up the particles with large local strain changes $\Delta\gamma_{i}=\gamma_{i}(B)-\gamma_{i}(A)>0.04$ between two points A and B at $T=0.2$ in (b), where $\bar{\gamma}$ is $0.016$ at B. We define the distribution function, $P_{\rm s}(v)=\sum_{i}{\delta(v-\Delta\gamma_{i})}/N$ for strain changes $\Delta\gamma_{i}=\gamma_{i}(t_{1})-\gamma_{i}(t_{0})$ between two times in (b). In (d), it is narrower for the initial interval ($\Delta\Phi:0\to 4$) than for the subsequent one ($\Delta\Phi:4\to 6$).

The shapes of our dipolar spheroids are centrosymmetric, leading to the electrostriction relation,

[TABLE]

at relatively high $T$. In Fig.14, Eq.(35) nicely holds with $C_{\rm es}\cong 0.012\sigma_{1}^{3}/\epsilon$ at $T=0.6$, while a closed loop appears at $T=0.2$. If we set $\epsilon/k_{\rm B}=100$ K and $\sigma_{1}=5~{}{\rm\AA}$, our $C_{\rm es}$ becomes 10 m${}^{4}/$C2. For ferroelectric polymers, Eq.(35) was found with a negative coefficient Zhang ; Z1 ; Furu1 ; Bl ($-13.5$ m${}^{4}/$C2 after electron irradiationZhang ). In contrast, the piezoelectric relation ($\bar{\gamma}\propto P_{z}$) holds for relaxors above the transitionpiezo .

VI ZFC$/$FC temperature changes

A large number of ZFC$/$FC experiments have been performed, where $T$ is varied at zero or fixed ordering field (electric fieldWest ; Uesu ; Z1 ; Ho1 , magnetic fieldF0 ; F1 ; F4 , and stressF5 ; F3 ). However, the physical pictures of these processes remain unclear. Here, we show relevance of collective, large-angle orientational changes in these cycles.

We followed cycles in Fig.15(a) setting $dT/dt=\pm 10^{-5}$ at fixed volume with $c=0.2$ and $\mu_{0}=0.8$. In ZFC-FH, (i) we cooled the system from a high-$T$ state to a low-$T$ state (point A) with $\Delta\Phi={P}_{z}=0$ and then (ii) heated it with $\Delta\Phi=0.5$ ($E_{\rm a}=0.024$) back to the initial $T$. Subsequently, in FC-FH, (iii) we cooled the system to point B with $\Delta\Phi=0.5$ and then (iv) heated it back with $\Delta\Phi=0.5$ fixed. We set $T=0.1$ at A and B.

In (a), we plot $\langle{n_{z}}\rangle={P}_{z}/\mu_{0}n_{1}$ vs $T$ on the two paths. The two heating curves meet at a freezing point C, where $T$ is given by $T_{\rm f}=0.21$ and the relaxation time $\tau_{1}$ in Eq.(27) is of order $10^{4}$. This $T_{\rm f}$ is very close to $T_{\rm m}$ at $\omega=2\pi\times 10^{-4}$ in Fig.6(a). Far below $T_{\rm f}$, the two curves are largely separated indicating marked nonergodicity, while they coincide for $T>T_{\rm f}$ in the ergodic regime.

In (b), we display $\langle{Q_{\ell}}\rangle=\sum_{i\in 1}Q_{\ell i}/N_{1}$ ($\ell=1,2$) in the same simulation run. From Eq.(24) they represent the average dipolar and quadrupolar orders. The difference of $\langle{Q_{1}}\rangle$ in the two cycles is at most $30\%$, while that of $\langle{Q_{2}}\rangle$ is only about $5\%$. Note that $Q_{2i}$ are rather insensitive to collective reorientations for most $i$ (see Sec.VB).

In (c) and (d), the dipoles with $Q_{1i}>0.7$ are depicted at A and B. These two patterns look similar, but some PNRs in the same locations in A and B have different polarization directions (for example, $[1\bar{1}\bar{1}]$ in A and $[\bar{1}11]$ in B). In the present example, the potential energy $U$ is $-5.86N\epsilon$ at A and $-5.88N\epsilon$ at B. Their difference ($=-0.02N\epsilon$) is small, but is still 5 times larger than $-E_{\rm a}M_{z}(=-0.004N\epsilon)$ at B (see $U_{\rm d}$ in Eq.(11)). Note that large potential barriers exist for reorientations of PNRs from the configurations at A to those at B. These barriers decrease with increasing $E_{\rm a}$, but its present size 0.024 is small. If a much larger $E_{\rm a}$ is applied at A, there can be a transition to a ferroelectric stateBl ; Bobnar .

In (e) and (f), we display the dipoles with large angle changes from A to C and from B to C. They satisfy $\Delta n_{zi}(A\to C)=n_{zi}(C)-n_{zi}(A)>0.8$ in (e) and $\Delta n_{zi}(B\to C)=n_{zi}(C)-n_{zi}(B)<-0.8$ in (f), where $n_{zi}$ are the $z$ component of {\mbox{\boldmathn}}_{i}. These large-angle changes are collective and heterogeneous. This should be a universal feature in glass coupled with a phase transition.

On the two FH paths, the potential barriers between the two states at the same $T$ remain very large for $T<T_{\rm f}$. They can be overcome by thermal activations at $T=T_{\rm f}$ (at C), where the reorientation rate of PNRs should be comparable to the inverse of the observation time $\tau_{\rm obs}$. Estimating the former as the inverse of $\tau_{1}$ in Eq.(27) and setting $\tau_{\rm obs}\sim T/(dT/dt)$, we obtain

[TABLE]

Indeed, we have $\tau_{1}=2.4\times 10^{4}\cong 1.1T/(dT/dt)$ at C. Here, $\tau_{1}$ at the freezing should decrease significantly for large $E_{\rm a}$ (not shown here). It follows that $T_{\rm f}$ at C decreases with increasing $\tau_{\rm obs}$. Note that this dependence is weak for long $\tau_{\rm obs}$ due to the abrupt $T$ dependence of $\tau_{1}$ at low $T$. It is well known that nucleation in a metastable state starts at an onset temperatureOnukibook , which is rather well defined for long $\tau_{\rm obs}$.

In Fig.15, we have used small $E_{\rm a}(=0.024$). However, cooling at a high electric field, Bobnar et al. Bobnar detected a field-induced ferroelectric transition. Such a phase transition can be predicted from a Ginzburg-Landau theory for $P_{z}$ including the field term $-E_{\rm a}P_{z}$ Bl ; Onukibook .

VII Antiferroelectric glass

So far we have treated ferroelectric glass. However, antiferroelectric order has been observed in mixtures containing cyanide units CN- such as KBr-KCN at low $T$ Mertz ; ori ; Binder ; anti-Loidl . It is also known that antiparallel alignment freezes at low $T$ in polar globular molecules such as cyanoadamantane anti ; anti0 containing CN or betaine phosphateAlbers containing H3PO4 due to their mutual steric hindrance. These systems should become antiferroelectric glass at low $T$ even without impurities. Here, we consider a mixture of dipoles and impurities introducing a short-range interaction favoring antiparallel ordering.

Supposing top-tail asymmetry of the dipoles, we replace the factor $A_{ij}$ ($i\in\alpha$ and $j\in\beta$) in Eq.(3) by

[TABLE]

The second term yields an exchange interaction between dipoles $i$ and $j$, where positive (negative) $J$ favors antiferroelectric (ferroelectric) ordering. We performed simulation for $J=0.1$ with $c=0.2$, $\mu_{0}=0.8$, and $\Delta\Phi=0$.

In Fig.16(a), we plot $\langle{Q_{1}}\rangle$ and $\langle{Q_{2}}\rangle$. Here, due to antiferroelectric ordering, $\langle{Q_{1}}\rangle$ remains very small at any $T$, but $\langle{Q_{2}}\rangle$ increases up to 0.42 with lowering $T$. Thus, the system exhibits quadrupolar order without dipolar order at zero applied electric field Binder ; ori . In more detail, we show the distribution functions $P_{\ell}(Q_{\ell})=\langle{\sum_{i\in 1}\delta(Q_{\ell}-Q_{\ell i})/N_{1}}\rangle$ at $T=0.05$. In (b), $P_{1}(Q_{1})$ is nearly symmetric (even) with respect to $Q_{1}\to-Q_{1}$ and $P_{2}(Q_{2})$ has a maximum at $Q_{2}\cong 1$ leading to $\langle{Q_{2}}\rangle\sim 0.4$.

Furthermore, in (c), a snapshot of the dipoles and the impurities is given, which looks very complicated. In (d), we show a typical antiferroelectric nanoregion in the middle of the cell, which are viewed from two directions. Any dipole $i$ in this region satisfies \mbox{\boldmathn}_{i}\cdot\mbox{\boldmathn}_{j}<-0.98 for some nearby $j$ with $r_{ij}<1.4$ within the same region. It is composed of 170 dipoles and surrounded by 130 impurities with no impurities in its interior. In (d) and (e), cross-sectional particle configurations are displayed at $z=H/2+2$ and $H/2+3$, respectively. We can see antiferroelectric ordering unambiguously for the dipoles parallel or antiparallel to the $z$ axis (perpendicular to $(111)$), while the orientations apparently look irregular for those perpendicular to $(\bar{1}11)$, $(1\bar{1}1)$, or $(11\bar{1})$.

VIII Summary and remarks

With molecular dynamics simulation, we have studied dipolar glass in mixtures of dipolar spheroids and apolar impurities in applied electric field. Properly calculating the electrostatics, we have visualized polar nanoregions (PNRs) and clarified their role in the dielectric response. We summarize our main results as follows.

(i) In Sec.II, we have introduced orientation-dependent Lennard-Jones potentials mimicking spheroidal repulsion. For its mild aspect ratio, the particles first form a fcc plastic crystal. Then, at lower $T$, the spheroids align along $\langle{111}\rangle$ resulting in rhombohedral structures. Assuming that each spheroid has a dipole parallel to its long axis, we have constructed an electrostatic energy $U_{\rm d}$ in Eq.(11), which accounts for the image dipoles and the applied field $E_{\rm a}$. In equilibrium, the differential susceptibility $\chi_{\rm dif}=dP_{z}/dE_{\rm a}$ is related to the polarization fluctuations as in Eq.(17).

(ii) In Sec.III, we have presented results on a structural phase transition in a one-component system of dipolar spheroids. It changes from a fcc crystal to a polycrystal with eight rhombohedral variants. This transition occurs in a narrow temperature range due to the finite size effect imposed by the metal walls.

(iii) In Sec.IV, we have examined diffuse ferroelectric transitions. The impurity distribution has been determined during crystallization, so marked impurity clustering has appeared. In our model, ferroelectric domains are broken up into smaller PNRs with increasing the impurity concentration $c$. For $c=0.2$, we have calculated the orientational time correlation function $C_{1}(t)$ in Fig.5 and the dielectric permittivity in Fig.6. The temperature of maximum of $\varepsilon^{\prime}$ is written as $T_{\rm m}(\omega)$. For very small $\omega$, the polarization fluctuations are enhanced for $T>T_{\rm m}$, but are composed of frozen PNRs and thermal fluctuations for $T<T_{\rm m}$. Individual PNRs have been visualized in Fig.10. The surface effects on the dipole orientations and the local electric fields have been examined in Sec.IVF and Appendix A.

(iv) In Sec.V, we have examined the polarization and the strain to cyclic applied electric field. At relatively high $T$, there is no hysteresis and an electrostriction relation holds. At low $T$, the polarization is on a hysteresis loop. In the cycle, collective large-angle changes are dominant where $\chi_{\rm dif}=dP_{z}/dE_{\rm a}$ is large.

(v) In Sec.VI, we have investigated the ZFC-FH and FC-FH thermal cycles in accord with the previous experiments. The frozen states at the lowest $T$ in the two cycles have been visualized in Fig.15. On the FH paths, heterogeneous collective reorientations have been found. These paths meet at a temperature $T_{\rm f}$, at which the reorientation rate ($\sim\tau_{1}^{-1})$ is of the same order as the ramping rate of the temperature ($\sim(dT/dt)/T$).

(vi) In Sec.VII, we have investigated antiferroelectric glass by introducing a short-ranged exchange interaction stemming from molecular shape asymmetry. We have visualized a typical antiferroelectric nanoregion.

(vii) In Appendix B, we have shown the method of calculating $\varepsilon^{\prime}$ and $\varepsilon^{\prime\prime}$ and found their algebraic behavior $(\propto\omega^{-\beta})$ at relatively large $\omega$ in the ergodic $T$ range.

Finally, we remark on future problems. (1) The isochoric specific heat $C_{V}$ can be calculated from the average energy. We found that it has a rounded peak in our mixture systems (not shown in this paper). This is consistent with the behavior of the isobaric specific heat $C_{p}$ in previous experiments ori ; Mertz ; Kawaji ; Tachibana . (2) There is a gradual crossover in the polarization fluctuations in the diffuse transition. For example, the PNRs have no clear boundaries at relatively high $T$, while sharp interfaces can appear at low $T$. It is of interest how the space correlations in the polarization and the particle displacements depend on $T$. (3) In real systems, impurities or mixed components have charges or dipoles. In solids, the polarization response can be large when ion displacements occur within unit cells as a phase transition. These features should be accounted for in future simulations. (4) Intriguing critical dynamics exists in the ergodic $T$ rangeKlee1 ; Bokov3 ; Bokov4 ; Cowley , as suggested by Eq.(31). The aging and memory effects at low $T$aging0 ; aging should also be studied in future (see the last paragraph of Sec.VB).

**Acknowledgments

**This work was supported by KAKENHI 15K05256, and KAKENHI 25000002. The numerical calculations were preformed on CRAY XC40 at YITP in Kyoto University and on SGI ICE XA/UV at ISSP in the University of Tokyo.

**Appendix A: Electrostatics of dipole systems

**

Here, we explain the electrostatics of dipoles between metal walls in applied field Takae1 ; Takae2 ; Hautman ; apply ; Klapp . The electric potential due to the image dipoles is equivalent to that due to the surface charge densities, written as $\sigma_{0}(x,y)$ at $z=0$ and $\sigma_{H}(x,y)$ at $z=H$. Without adsorption and ionization on the surfaces, the dipole centers are somewhat away from the walls (see the comment below Eq.(6)). Then,

[TABLE]

where $E_{z}=-\partial\Phi/\partial z$. We consider the 2D Fourier expansions of $\sigma_{\lambda}$. For $\lambda=0$ and $H$ they are

[TABLE]

where {\mbox{\boldmathr}}_{\perp}=(x,y), and {\mbox{\boldmathk}}=(2\pi/L)(n_{x},n_{y})\neq(0,0) with $n_{x}$ and $n_{y}$ being integers. The first term is the mean surface charge density ${\bar{\sigma}}_{\lambda}=\int_{0}^{L}dx\int_{0}^{L}dy~{}\sigma_{\lambda}(x,y)/L^{2}$. From Eq.(10) we can express the Fourier components \sigma_{\lambda{\mbox{\boldmathk}}} asapply

[TABLE]

where \nabla_{j}=\partial/\partial{\mbox{\boldmathr}}_{j}, $G_{k}^{0}(z)={\sinh(k(H-z))}/\sinh(kH)$, and $G_{k}^{H}(z)={\sinh(kz)}/\sinh(kH)$ with k=|{\mbox{\boldmathk}}|.

For dipolar systems, the Poisson equation is written as

[TABLE]

Integration of Eq.(A4) in the cell yields ${\bar{\sigma}}_{0}+{\bar{\sigma}}_{H}=0$. We also multiply Eq.(A4) by $z$ and integrate it in the cell. Using the total polarization $M_{z}$ we findHautman ; Takae1

[TABLE]

without surface adsorption and ionization. The fluctuations of ${\bar{\sigma}}_{0}$ and $M_{z}/V$ thus coincide at fixed $E_{\rm a}$ Hautman ; apply ; Takae1 ; Takae2 .

The mean surface charge densities produce the potential $-4\pi{\bar{\sigma}}_{0}z$ in the cell, so $\Phi$ consists of three parts as

[TABLE]

which is equivalent to Eq.(10). The first term $\Phi_{\rm d}$ arises from the dipoles in the cell. Imposing the lateral periodic boundary condition, we express it as

[TABLE]

where {\mbox{\boldmathg}}({\mbox{\boldmathr}})=r^{-3}{\mbox{\boldmathr}} and {\mbox{\boldmathm}}_{\perp}=(m_{x},m_{y},0) with $m_{x}$ and $m_{y}$ being integers. The third term $\phi_{\rm s}$ in Eq.(A6) arises from the charge density deviations $\delta\sigma_{\lambda}(x,y)=\sigma_{\lambda}-{\bar{\sigma}}_{\lambda}$. In terms of \sigma_{\lambda{\mbox{\boldmathk}}} in Eq.(A2), $\phi_{\rm s}$ is expressed as

[TABLE]

Now the local electric field {\mbox{\boldmathE}}_{i} is written as

[TABLE]

The first term arises from the other dipoles in the cell:

[TABLE]

The second term is due to the surface charges:

[TABLE]

where the first term is homogeneous and {\mbox{\boldmathE}}_{\rm s}({\mbox{\boldmathr}})=-\nabla\phi_{\rm s}({\mbox{\boldmathr}}) is due to $\delta{\sigma_{\lambda}}={\sigma_{\lambda}}-{{\bar{\sigma}}_{\lambda}}$. The dipoles next to the walls are parallel or antiparallel to the $z$ axis due to {\mbox{\boldmathE}}_{\rm s} even for $\Delta\Phi=0$ (see the snapshots in this paper)apply . However, as in Fig.17(b), {\mbox{\boldmathE}}_{\rm s} is negligibly small (even in ferroelectric states) if the distances from the walls exceed the typical domain size. This is due to the factors $\exp({-kz})$ and $\exp(-k(H-z))$ in Eq.(A8).

In Fig.17, we show (a)$\sigma_{0}(x,y)$, (b)$\overline{E_{\rm s}^{2}}(z)$, (c)$E_{\rm s}(x,y,1)$, and (d)$E_{\rm s}(x,y,3)$ in ferroelectric glass of our system, where we set E_{\rm s}(x,y,z)=|{\mbox{\boldmathE}}_{\rm s}|=|\nabla\phi_{\rm s}| and

[TABLE]

Here, $\sigma_{0}(x,y)$ in (a) and $E_{\rm s}(x,y,1)$ in (c) consist of microscopic and mesoscopic fluctuations. The latter arise from the PNRs near the surface from Fig.11, being apparent in (d). For $z$ longer than the PNR length, $E_{\rm s}(x,y,z)$ decays to zero in (b). Thus, {\mbox{\boldmathE}}_{i}^{\rm sur}\to-4\pi{\bar{\sigma}}_{0}{\mbox{\boldmathe}}_{z} far from the walls, which was previously found for liquid waterTakae1 ; Takae2 .

**Appendix B: Linear response to oscillating field and frequency-dependent susceptibilities

**

We applied a small sinusoidal electric field of the form $E_{\rm a}(t)=E_{\rm a}^{1}\sin(\omega t)$ with $E_{\rm a}^{1}=0.047$. We calculated the polarization response $\delta P_{z}=\delta M_{z}/V$ to this perturbation over 10 periods. After a few periods, it is expressed as

[TABLE]

where $\chi^{\prime}$ and $\chi^{\prime\prime}$ are the frequency-dependent susceptibilities. Then, $\varepsilon^{\prime}$ and $\varepsilon^{\prime\prime}$ in Fig.6 are defined by

[TABLE]

The Hamiltonian ${\cal H}$ increases as ${\overline{d{\cal H}/dt}}=V\omega\chi^{\prime\prime}|E_{\rm a}^{1}|^{2}/2$ in time (see below Eq.(20)), where the time average is taken in one period. From Eq.(A5) the mean surface charge density at the bottom wall is written as

[TABLE]

which oscillates as $\sin(\omega t-\delta_{\rm p})$ with $\tan\delta_{\rm p}=\varepsilon^{\prime\prime}/\varepsilon^{\prime}$. In Fig.6, we give the resultant $\varepsilon^{\prime}$ and $\varepsilon^{\prime\prime}/\varepsilon^{\prime}$ in a wide $T$ range including the nonergodic range.

On the other hand, around equilibrium, we can use the linear response theory Kubo for the Hamiltonian of the form (15). Within this scheme, the dielectric response can be expressed in terms of the time-correlation function for the deviation $\delta M_{z}(t)=M_{z}(t)-{\langle{M_{z}}\rangle}_{\rm e}$:

[TABLE]

where ${\langle{\cdots}\rangle}_{\rm e}$ represents the equilibrium average and $\chi={\langle{(\delta M_{z})^{2}}\rangle}_{\rm e}/Vk_{\rm B}T=(\varepsilon-1)/4\pi$ at $E_{\rm a}=0$ (see Eq.(17)). Using $G_{z}(t)$ we obtain the linear response relations,

[TABLE]

The complex susceptibility $\chi^{\prime}-i\chi^{\prime\prime}$ can be expressed as

[TABLE]

in terms of the time derivative ${\dot{G}}_{z}(t)={dG_{z}(t)}/{dt}$.

In Fig.18(a), we show our numerical results of $G_{z}(t)$ in Eq.(B4) at $T=0.5,0.6,0.7$, and 1.0 for $c=0.2$, $\Delta\Phi=0$, and $\mu_{0}=1.6$, where the data at long times are inaccurate, however. From $G_{z}(\tau_{G})=e^{-1}$, we define the relaxation time $\tau_{G}$, which is somewhat shorter than $\tau_{1}$ in Fig.5. In fact, we obtain $(\tau_{G},\tau_{1})=(230,1300),(36,60)$, $(11,17)$, and $(1.8,2.5)$ for $T=0.5,0.6$, 0.7, and 1.0, respectively. We may introduce another time by $\tau^{\prime}_{G}=\lim_{\omega\to 0}\chi^{\prime\prime}/\omega\chi=\int_{0}^{\infty}dtG_{z}(t)$, but we confirm $\tau_{G}^{\prime}\sim\tau_{G}$.

In (b), the initial decay of $G_{z}(t)$ is well fitted to

[TABLE]

where $\beta\cong 0.57$ and $A_{\rm p}\cong 0.63$. Then, ${\dot{G}}_{z}(t)\propto-t^{\beta-1}$ for $t\lesssim\tau_{G}$. If this is substituted into Eq.(B7), we find

[TABLE]

where $B_{\rm p}=\beta\Gamma(\beta)A_{\rm p}\chi\sim\chi$. The algebraic form (B9) with $0<\beta<1$ has been observed in many systems including relaxors and mixed crystals Bokov3 ; Klee ; ori .

We calculated the ratios $\chi^{\prime}(\omega)/\chi$ and $\chi^{\prime\prime}(\omega)/\chi$ from Eqs.(B5) and (B6) using $G_{z}(t)$ in (a). In (c), they are plotted at $T=0.6$ together with those from the data in Fig.6, where the latter are from Eq.(B1). The points from these two sets fairly agree for any $\omega$. In (d), we also plot these ratios vs $\omega\tau_{G}$ at three temperatures using the results in Fig.6. The behaviors in the region $\omega\tau_{G}\gtrsim 1$ in (c) and (d) support the algebraic form (B9).