Electromechanical properties of ferroelectric polymers: Finsler geometry modeling and Monte Carlo study

Monte Carlo simulations are performed to study the electromechanical properties of polyvinylidene difluoride (PVDF) because of its unusual deformation in response to an applied electric field. PVDF is a ferroelectric polymer characterized by a negative strain along the applied electric field direction. However, the mechanism of this electromechanical response remains unclear due to the complex nature arising from the hierarchical structure across length scales. In this paper by employing the Finsler geometry (FG) model as a new approach to this problem, we show that the deformations observed in the 3D tetrahedral model are almost identical to those of real PVDF. More precisely, the simulated mechanical deformation and polarization, which are the responses to the electric field, are almost exactly the same as those in real experiments.

Introduction. -Ferroelectric polymers, such as polyvinylidene difluoride (PVDF), are known to have strong coupling between their mechanical and electrical properties. Indeed, PVDF and its copolymers are commonly used for designing actuators because of their large deformations under external electric fields [1][2][3]. In contrast to typical classical ferroelectrics, PVDF has an unusual negative longitudinal piezoelectric coefficient, i.e., the bulk thickness becomes thinner in the applied electric field direction. The mechanical strain is well known to be proportional to the square of the polarization [4]. This implies that the field-induced deformation in PVDF mainly comes from electrostriction. However, the relation between the strain and polarization under low electric fields is reported to be linear [4]. This complex nature comes from the fact that both piezoelectricity and electrostriction are associated with scales ranging from the monomer size to the polymeric domain size, and this complex multiscale behavior can be summarized by the "elec-tromechanical properties" of PVDF.
Broadhurst et al. [5] demonstrated that the origin of piezoelectricity can be explained by the mechanism called the dimensional effect. This mechanism assumes that the dipoles are rigid and attributes polarization to macroscopic shape deformation. Katsouras et al. [6] revealed that PVDF piezoelectricity comes from two microscopic mechanisms: changes in the lattice constant and coupling of crystalline and amorphous parts.
A number of studies based on first principles and molecular dynamics simulations have been devoted to PVDFbased polymers [7][8][9][10]. However, purely atomistic approaches are not always efficient for studying polymers because polymer chains consist of thousands of monomers, which is beyond the limit of such calculations. In contrast, well-developed phenomenological models for ferroelectrics exist that can reproduce polarization-electric field (PE) curves [11]. Several works in which molecular dynamics and phenomenological approaches are combined have also p-1 been presented [12]. However, despite such a long history of studies, the molecular mechanism of the piezoelectric effect and the electrostriction in PVDF still remains unclear.
In the Finsler geometry (FG) model, materials are assumed to be a Finsler space rather than in a Finsler space [13]. In this approach, the anisotropy of the mechanical, optical or other properties is represented by a Finsler metric, which effectively deforms discrete differential operators. As a consequence, the corresponding interaction between molecules, described by discrete differentials, becomes direction dependent and hence anisotropic. The FG modeling technique has been successfully applied for evaluating the deformation of materials with mechanical anisotropy, such as rubbers and soft biological materials [14]. Moreover, the shape transformation of liquid crystal elastomers under external electric fields has been successfully studied by this technique [15]. In this Letter, we further extend the FG model to describe the unusual piezoelectric effect in ferroelectric polymers, where the Finsler metric is defined by using an internal degree of freedom σ corresponding to the polarization [16].
The model. -We begin with the continuous form of the Hamiltonian for the tensile energy where r(∈ R 3 ) is a position vector of the material, and x a (a = 1, 2, 3) are the local coordinates. The symbol g ab is the inverse of Finsler metric g ab , which will be given below, and g is its determinant. If g ab is the Euclidean metric, then the expression (1) becomes the classical Gaussian bond potential, which is a 3D extension of the linear chain model [17].
To define the Finsler metric, let us consider a 3D thin cylinder discretized by tetrahedrons using Voronoi tessellation ( Fig. 1(a)). The ratio of the cylinder height and its diameter is 0.125 (see [15] for more details). A variable σ i (∈ S 2 : unit sphere) is introduced for the polarization at vertex i of a tetrahedron ( Fig. 1(b)). Using this σ i , we define the Finsler length along bond ij such that where t ij is the unit tangential vector from vertex i to vertex j. The expression in Eq. (2) is new and different from that used in [15] and defined such that tetrahedrons shrink along the direction of σ. The parameter v 0 plays a role in the strength of the anisotropy in the mechanical properties, and v 0 = 0.1 is assumed in this study. At vertex 1 of tetrahedron 1234 ( Fig. 1(b)), the metric is defined by Using the equivalence of expressions with respect to sub- , · · · and replacing differentials with differences, we obtain whereN is the mean total number of tetrahedrons sharing bond ij and is given byN ≃ 4.84 for the lattice used here. The full Hamiltonian of the model including some additional terms is given by We should note that the coefficient γ of S 1 is fixed to γ = 1 for simplicity. S 2 plays a role in the strength against shear and bending deformations, and φ i denotes the internal angle of the triangle at vertex i ( Fig. 1(b)). The parameter κ is a stiffness and plays a role in the strength of the resistance against all deformations of the tetrahedron except the similarity one. We should note that, as we will see later, κ effectively determines the strength of the electromechanical coupling, i.e., it is macroscopically reflected in the slope of the strain-polarization curve.
S 0 represents an interaction between two nearest neighbors σ i and σ j . Such a Heisenberg spin Hamiltonian is widely used in simulations for ferroelectric domain structures [18][19][20]. In those models, the energy σ i · σ j = −S 0 /N without an external electric field increases with increasing λ. This is almost the same as in our FG model except that σ interacts with the lattice deformation via v ij p-2 in Eq. (3). Thus, the remnant polarization is controlled by λ.
To describe the interaction between the polarization and external electric field E, we introduce S 3 and S 4 . Note that the coefficient for S 3 is assumed to be 1 for simplicity. S 3 is a classical PE field interaction, and S 4 is its quadratic extension. As a consequence, the shape of the PE curve is expected to be dependent on the parameter α.
Indeed, for small α, the PE hysteresis loop has a squarelike shape, which is typical of crystalline materials. The loop becomes more rounded with increasing α. Notably, λ and α determine the polarization strength of the model, while the electromechanical coupling can be controlled by varying the parameters κ and v 0 .
Simulation results. -The Monte Carlo (MC) simulation technique is used to study the behavior of the 3D cylinder under external field E along the z axis ( Fig. 1(a)). The vertex position r and polarization σ are updated by the Metropolis algorithm [21]. To suppress meaningless configurations, we introduce several constraints on r. The volume of a tetrahedron is prohibited from being negative and from exceeding tens of its initial mean value. The squared bond lengths are constrained to be less than 10ℓ 2 0 , where ℓ 0 is the mean initial bond length. To keep the thin cylinder horizontal, a hard-wall potential is introduced so that the z component of r is inside (−H 0 , H 0 ), where H 0 is the mean initial height of the cylinder in the equilibrium configuration for E = 0 and is evaluated in the test simulations.
In this Letter, we compare our simulation results with experimental data of uniaxially drawn PVDF measured at a low electric field frequency of 0.1 Hz [4]. In those experimental data, the polarization and strain vs. electric field always form hysteresis loops. However, the standard MC technique only allows us to study equilibrium properties. Therefore, we simulate only parts of the return loop where the polarization aligns along the electric field and switching is not expected, as in the equilibrium configuration. Indeed, the data acquisitions were done at frequencies less than 1 kHz [22], and hence, 0.1 Hz is sufficiently slow for equilibration.
In the simulations, we fix λ = 0.345 to make the system be in the ferroelectric phase. The reason for the choice of λ = 0.345 is as follows: According to the test simulations, the polarization σ z starts to saturate at σ max z = 0.94 with further increases in the electric field. Additionally, in the experimental data, the maximal experimental polarization is P max exp = 0.097 C/m 2 = 97 mC/m 2 [4]. Hence, σ z and P exp are related by σ z = ( σ max z /P max exp )P exp ≈ 9.7P exp . This relation can also be used for the remnant polarization P r = 54.8 mC/m 2 , which is the experimental polarization as E → 0 [4]. Thus, we obtain the simulated remnant polarization at E = 0 such that σ r z = 0.0548 × 9.7 = 0.53, which leads us to λ = 0.345 in the test simulations. Thus, using this factor P r / σ r z = 54.8/0.53 ≈ 103.4, we define P by P = 103.4 σ z to com-pare σ z with P exp .
Next, we consider the relation between the external electric fields E and E exp . First, we have [15] ǫ 0 ∆ǫE 2 exp = αE 2 k B T a 3 (6) for the electrostriction energy, where ǫ 0 = 8.85×10 −12 F/m is the vacuum permittivity, ∆ǫ = 0.25 is the dielectric anisotropy taken from [23], and a is the lattice spacing. The parameter α is fixed relatively large (α = 300) because PVDF is a semicrystalline polymer, in which some parts of the material are in the amorphous state. This is also the reason why the shape of the PE loop becomes more rounded compared to crystalline ferroelectrics. From the comparison of the dipole-field interactions, we also have where f is understood to be the number of monomers associated with one vertex, N = 10346 is the total number of vertices, and V = 836 is the cylinder volume (in units of a 3 ) at E = 0.15. From Eqs. (6) and (7), we have E = ǫ 0 ∆ǫ(N σ z /αV P exp )f E exp ≃ 7.96 × 10 −10 E exp using the abovementioned numbers if f = 900, where E exp is given in units of MV/m. At the same time, we have a ≈ 7.01×10 −9 m from Eq. (6) using k B T = 4 × 10 −21 . This a is much larger than the van der Waals distance 10 −10 m, i.e., a lies in a reasonable range. The number f = 900 is also reasonable. This f is slightly larger than the estimated real density (N A ρV a 3 )/(M N ) ≈ 570, where ρ = 1.97 g/(cm) 3 is the density of β-PVDF, N A is Avogadro's number, and M (=64 g/mol) is the molar mass of the PVDF monomer. The reason for this deviation is that our model is coarse-grained, and each vertex corresponds to a lump of monomers; consequently, the monomer density can be different from the actual density. Moreover, if the coefficient γ of γS 1 in Eq. (5) is taken to be slightly larger than 1, then the volume will be slightly smaller than the current V , and consequently, f becomes smaller and is hence controllable.
First, snapshots of the PVDF model are shown in Figs.  2(a),(b), where the external electric field is E = 0 and E = 0.15, which corresponds to E exp ≃ 189 MV/m. The small red cylinders on the surface denote the variable σ. σ aligns along the z axis even when E = 0 because λ is fixed to λ = 0.345, corresponding to the ferroelectric behavior of the β phase of PVDF.
Based on Fig. 3(a), the simulation results of P vs. E are in good agreement with P exp vs. E exp . The influence of other simulation parameters κ and v 0 , which are fixed to κ = 1.8 and v 0 = 0.1, on the PE curve is almost negligible.
In Fig. 3(b), the height strain ε H vs. electric field is presented. The simulated strain is obtained by measuring the mean height of the cylinder. The mean height is calculated from the vertices on the upper and lower surfaces at the central part inside a circle of radius R 0 = 10ℓ 0 , сϬ сϬ͘ϭϱ where ℓ 0 is the mean bond length. We should note that the experimental loop is of the butterfly shape, and the polymer shrinks (elongates) along the electric field when the field direction and polarization are parallel (antiparallel) to each other. The simulation results show that the height of the cylinder shrinks with increasing electric field, in agreement with the experimentally observed data.
The dependence of the PVDF strain on the square of the polarization is linear and has almost no hysteresis [4]. This implies that the strain is directly related to the polarization rather than to the electric field. Therefore, we expect that the electrostriction constant, i.e., the slope of the strain-polarization curve, can be controlled by κ. The solid line in Fig. 4(a), denoted "Theory", is drawn by using P 2 (S) = P 2 (0) + (1/κ 33 )S with the electrostriction constant κ 33 = −2.4 and the experimental remnant polarization P (0), where S denotes the strain. The constant |κ 33 | decreases with increasing κ. The experimental sample has the constant κ 33 = −2.4, and this value is achieved in our model with κ = 1.8.
The elongation of the cylinder diameter D is also calculated by the formula D = 2 V /(πH 0 ) with the initial height H 0 at E = 0. Using this D, we obtain the diameter strain ε D , which is considered to be the "nominal strain" because H 0 is used for calculating D; hence, we estimate the Poisson ratio ν =−ε D /ε H (Fig. 4(b)). ν tends to be close to 0.3 with increasing strain inside the experimental electric field range. This value is comparable to the experimental value of 0.33 for β-PVDF obtained by mechanical testing [24].
Thus, we conclude that the FG model successfully reproduces the electromechanical properties of β-PVDF. We must emphasize that these PE, strain-electric (SE) and Poisson ratio data are obtained with a single set of parameters.
Concluding remarks. -We have shown that the FG model is applicable to the study of the electromechanical properties of ferroelectric polymers. The Monte Carlo simulation results are in good agreement with the experimental data of β-PVDF. Both PE field and SE field curves are reproduced with a single set of simulation parameters. Moreover, the simulated Poisson ratio is also reasonable. The technique used for PVDF can also be applicable to ferroelectric ceramics such as BaTiO 3 and to magnetostriction of ferromagnetic materials, which will be studied in the future. * * *