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

Ferroelectric polymers such as polyvinylidene difluoride (PVDF) exhibit strongly coupled mechanical and electrical properties. Because of their large deformations under external electric fields, PVDF and its copolymers are commonly used to design actuators [1, 2, 3]. In contrast to typical classical ferroelectrics, PVDF exhibits an unusual negative longitudinal piezoelectric coefficient, i.e., the bulk thickness decreases in the direction of the applied electric field. Moreover, the mechanical strain is proportional to the square of the polarization [4], which implies that the field-induced deformation in the PVDF is a result of the electrostriction. However, the strain and polarization are linearly correlated under low intensity electric fields [4]. This complex behavior occurs because both piezoelectricity and electrostriction are associated with scales ranging from the monomer to polymeric domain sizes, and this complex multiscale behavior can be reflected through the “electromechanical properties” of PVDF.

Broadhurst et al. [5] demonstrated that the piezoelectricity can be attributed to a mechanism known as the dimensional effect. According to this mechanism, the dipoles are rigid, and the polarization behavior can be attributed to the macroscopic shape deformation. Katsouras et al. [6] indicated that the piezoelectric nature of PVDF is a result of two microscopic mechanisms: changes in the lattice constant and coupling of the crystalline and amorphous parts.

Several studies based on first principles and molecular dynamics simulations have been devoted to PVDF-based polymers [7, 8, 9, 10]. However, purely atomistic approaches cannot be used to effectively examine polymers because polymer chains consist of thousands of monomers that cannot be processed through these calculations. Nevertheless, well-developed phenomenological models for ferroelectrics exist, which can reproduce the polarization-electric field (PE) curves [11]. In addition, several researchers combined molecular dynamics and phenomenological approaches [12]. Despite the notable research conducted in this domain, the molecular mechanisms of the piezoelectric effect and electrostriction in PVDF remain unclear.

In Finsler geometry (FG) modeling, the anisotropy of mechanical, optical or other properties is represented through a Finsler metric that effectively deforms the discrete elastic energy [13]. Consequently, the corresponding interaction among the molecules is direction dependent and hence anisotropic. The FG modeling technique has been successfully applied to evaluate the deformation of materials with mechanical anisotropy, such as rubbers and soft biological materials [14] as well as the shape transformation of liquid crystal elastomers under external electric fields [15]. In this letter, we extend the FG model to describe the unusual piezoelectric effect pertaining to ferroelectric polymers by defining the Finsler metric using an internal degree of freedom $\sigma$ corresponding to the polarization [16].

Consider the continuous form of the Hamiltonian for the tensile energy,

where $\vec{r}(\in{\bf R}^{3})$ is a position vector of the material and $x^{a}(a\!=\!1,2,3)$ denotes the local coordinates. The symbol $g^{ab}$ indicates the inverse of the Finsler metric $g_{ab}$, as described in the following text, and $g$ is the corresponding determinant. If $g_{ab}$ is the Euclidean metric, the expression (1) represents the classical Gaussian bond potential, which is a 3D extension of the linear chain model [17].

To define the Finsler metric, we consider a 3D thin cylinder discretized by tetrahedrons by using Voronoi tessellation (Fig. 1(a)). The ratio of the cylinder height to its diameter is 0.125, and the lattice is the same as that used in Ref. [15]. A variable ${\vec{\sigma}_{i}}~{}(\in S^{2}:{\rm unit\;sphere})$ is introduced for the polarization at vertex $i$ of a tetrahedron (Fig. 1(b)). In contrast to the case of a liquid crystal elastomer, as described in Ref. [15], ${\vec{\sigma}_{i}}$ is assumed to be a polar variable. Using ${\vec{\sigma_{i}}}$, we define the unit Finsler length along bond $ij$ such that

where $\vec{t}_{ij}$ is the unit tangential vector from vertex $i$ to vertex $j$. The expression presented as Eq. (2) is new and different from that used in [15]. In accordance with this new definition, the tetrahedrons shrink along the direction of ${\vec{\sigma}}$. The parameter $v_{0}$ serves to strengthen/weaken the anisotropy in the mechanical properties, and it is assumed that $v_{0}\!=\!0.1$. At vertex 1 of tetrahedron 1234 (Fig. 1(b)), the metric is defined as

Using the equivalence of the expressions with respect to the substitutions of indices $(1\rightarrow 2,2\rightarrow 3,3\rightarrow 4,4\rightarrow 1)$, $(1\rightarrow 3,2\rightarrow 4,3\rightarrow 1,4\rightarrow 2)$, $\cdots$ and by replacing the differentials with differences, we obtain

where $\bar{N}$ is the mean total number of tetrahedrons sharing bond $ij$ and is given by $\bar{N}\!\simeq\!4.84$ for the considered lattice.

We intuitively demonstrate the deformation mechanism of the shape of the liquid crystal elastomer by the external electric field. Note that $\Gamma_{ij}$ in $S_{1}$ functions as the local tension coefficient of the bond $ij$. If the variable $\sigma$ is aligned through an external electric field, owing to the interaction implemented in $v_{ij}$ in Eq. (2), almost all $\Gamma_{ij}$ along the direction of the electric field become larger than those along the perpendicular direction. Consequently, the corresponding bond lengths $\ell_{ij}$ along the field direction decrease. These aspects pertain to an intuitive explanation of the effect of the interaction of the polarization vector $\sigma$ and polymer position ${\vec{r}}$ implemented in $v_{ij}$.

The full Hamiltonian of the model including certain additional terms can be expressed as

The unit of energy is $k_{B}T$, which is fixed as $k_{B}T\!=\!1$ in the simulations, where $k_{B}$ and $T$ denote the Boltzmann constant and temperature, respectively. Notably, the simulation unit is defined by the relation $k_{B}T\!=\!1$ for the energy, and the relation for the lattice spacing $a\!=\!1$ pertains to the length. The lattice spacing $a$ is introduced such that all the quantities with the unit of length are multiplied by $a$. Moreover, the coefficient $\gamma$ of $S_{1}$ is fixed as $\gamma\!=\!1$ for simplicity. This condition can be ensured by rescaling $T$ in the Boltzmann factor $\exp(-S/k_{B}T$) such that $\gamma/T\!=\!1/T^{\prime}$. In accordance with the new $T^{\prime}$, all other coefficients are rescaled. The symbol $T$ denotes the temperature instead of $T^{\prime}$, and the expression for $S$ in Eq. (12) can be obtained. $S_{2}$ corresponds to the strength against the shear and bending deformations, and $\phi_{i}$ denotes the internal angle of the triangle at vertex $i$ (Fig. 1(b)). The parameter $\kappa$ denotes the stiffness, and it influences the resistance against all deformations of the tetrahedron except the similarity deformation. As described in the subsequent sections, $\kappa$ effectively determines the strength of the electromechanical coupling; i.e., this parameter is macroscopically reflected in the slope of the strain-polarization curve.

$S_{0}$ represents the interaction between two nearest neighbors $\sigma_{i}$ and $\sigma_{j}$. Such a Heisenberg spin Hamiltonian is widely used in simulations for ferroelectric domain structures [18, 19, 20]. In these models, the energy $\langle\vec{\sigma}_{i}\cdot\vec{\sigma}_{j}\rangle\!=\!-S_{0}/N$ without an external electric field increases with increasing $\lambda$. This phenomenon is the same as that considered in our FG model except $\sigma$ interacts with the lattice deformation through $v_{ij}$ in Eq. (3). Thus, the remnant polarization is controlled by $\lambda$.

To describe the interaction between the polarization and external electric field $\vec{E}$, we introduce $S_{3}$ and $S_{4}$. Note that the coefficient for $S_{3}$ is assumed to be 1 for simplicity. This condition can be ensured by rescaling the electric field from $E$ to $E^{\prime}$ such that $aS_{3}(E)\!+\!bS_{4}(E)\!=\!S_{3}(E^{\prime})\!+\!\alpha S_{4}(E^{\prime})$ for any combination of $a(\not=\!0)$ and $b$. Using the relations $E^{\prime}\!=\!aE$ and $\alpha\!=\!b/a^{2}$, we can obtain the expression $S_{3}\!+\!\alpha S_{4}$ in Eq. (12). The symbol $E$ denotes the electric field instead of $E^{\prime}$. $S_{3}$ represents a classical PE field interaction, and $S_{4}$ is the corresponding quadratic extension. In this case, the shape of the PE curve is expected to be dependent on the parameter $\alpha$. For small $\alpha$ values, the PE hysteresis loop is squarelike, which is typical of crystalline materials. The loop becomes more rounded with increasing $\alpha$. Notably, $\lambda$ and $\alpha$ determine the polarization strength of the model, and the electromechanical coupling can be controlled by varying the parameters $\kappa$ and $v_{0}$.

The partition function $Z$ is defined as

where $\sum_{\sigma}$ denotes the sum over all possible configurations of $\sigma$, and $\int\prod_{i=1}^{N-1}d{\vec{r}}_{i}$ denotes $3(N-1)$-dimensional integration, where the center of mass is fixed at the origin of ${\bf R}^{3}$.

The Monte Carlo (MC) simulation technique is used to study the behavior of a 3D cylinder under an external field $\vec{E}$ along the $z$ axis (Fig. 1(a)). The vertex position ${\bf r}$ and polarization $\vec{\sigma}$ are updated with the Metropolis algorithm [21]. To suppress meaningless configurations, we apply several constraints on ${\bf r}$. Specifically, the volume of a tetrahedron cannot be negative and exceed tens of the initial mean value. The squared bond lengths must be less than $10\ell_{0}^{2}$, where $\ell_{0}$ is the mean initial bond length. To ensure that the thin cylinder is horizontal, a hard-wall potential is introduced such that the $z$ component of ${\bf r}$ lies in $(-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 through test simulations. These constraints do not influence the equilibrium configurations, and this aspect can be verified by the relation $S_{1}/N\!=\!3(N-1)/(2N)\!\simeq\!1.5$, the plot of which is presented below. This relation can be attributed to the scale-invariant property of the partition function [15] and is used to verify whether the simulation program and simulations are sufficiently correct.

In this letter, we compare our simulation results with the experimental data of uniaxially drawn PVDF measured at a low electric field frequency of 0.1 Hz [4]. In the experimental data, the polarization and strain vs. electric field form hysteresis loops. However, the standard MC technique allows the examination of only the equilibrium properties. Therefore, we simulate only parts of the return loop in which the polarization aligns along the electric field and switching is not expected, as in the equilibrium configuration. The data are acquired at frequencies less than 1 kHz [22], and hence, a value of 0.1 Hz is considered to be sufficiently low for equilibration.

In the simulations, we set $\lambda\!=\!0.345$ to ensure that the system is in the ferroelectric phase. $\lambda\!=\!0.345$ is selected owing to the following: According to the test simulations, the polarization $\langle\sigma_{z}\rangle$ starts to saturate at $\langle\sigma_{z}^{\rm max}\rangle\!=\!0.94$ with a further increase in the electric field intensity. Additionally, among the experimental data, the maximal experimental polarization is $P_{\rm exp}^{\rm max}\!=\!0.097\text{ C}/\text{m}^{2}\!=\!97\text{ mC}/\text{m}^{2}$ [4]. Hence, $\langle\sigma_{z}\rangle$ and $P_{\rm exp}$ are related as $\langle\sigma_{z}\rangle\!=\!(\langle\sigma_{z}^{\rm max}\rangle/P_{\rm exp}^{\rm max})P_{\rm exp}\!\approx\!9.7P_{\rm exp}$. This relation can also be used to determine the remnant polarization $P^{\rm r}\!=\!54.8$ $\text{ mC}/\text{m}^{2}$, which corresponds to the experimental polarization as $E\!\to\!0$ [4]. Thus, we obtain the simulated remnant polarization at $E\!=\!0$ such that $\langle\sigma_{z}^{\rm r}\rangle\!=\!0.0548\!\times\!9.7\!=\!0.53$, and thus, $\lambda\!=\!0.345$ in the test simulations. Using the factor $P^{\rm r}/\langle\sigma_{z}^{\rm r}\rangle\!=\!54.8/0.53\!\approx\!103.4$, we define $P$ by $P\!=\!103.4\langle\sigma_{z}\rangle$ to compare $\langle\sigma_{z}\rangle$ with $P_{\rm exp}$.

Next, we consider the relation between the external electric fields $E$ and $E_{\rm exp}$. First, [15]

for the electrostriction energy, where $\epsilon_{0}\!=\!8.85\!\times\!10^{-12}$ F/m is the vacuum permittivity, $\Delta\epsilon\!=\!0.25$ is the dielectric anisotropy, referenced from [23], and $a$ is the lattice spacing. The parameter $\alpha$ is set to be relatively large ($\alpha\!=\!300$) because PVDF is a semicrystalline polymer, and thus, certain parts of the material are in the amorphous state. Moreover, owing to this reason, the PE loop becomes more rounded compared to that for the crystalline ferroelectrics. According to the comparison of the dipole-field interactions,

where $f$ is 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$. According to Eqs. (14) and (15), $E\!=\!\epsilon_{0}\Delta\epsilon(N\langle\sigma_{z}\rangle/\alpha VP_{\rm exp})fE_{\rm exp}\!\simeq\!7.96\!\times\!10^{-10}E_{\rm exp}$, using the abovementioned values, with $f\!=\!900$, where $E_{\rm exp}$ is in units of MV/m. Furthermore, $a\!\approx\!7.01\!\times\!10^{-9}$ m according to Eq. (14), considering $k_{B}T\!=\!4\!\times\!10^{-21}$. This $a$ is considerably larger than the van der Waals distance $10^{-10}$ m; i.e., $a$ lies in a reasonable range.

Moreover, the value $f\!=\!900$ is reasonable. This $f$ is slightly larger than the estimated real density $(N_{A}\rho Va^{3})/(MN)\!\approx\!570$, where $\rho\!=\!1.97\;{\rm g/(cm)^{3}}$ is the density of $\beta$-PVDF, $N_{A}$ is Avogadro’s number, and $M$(=64 g/mol) is the molar mass of the PVDF monomer. This deviation occurs because the model is coarse grained, and each vertex corresponds to a lump of monomers; consequently, the monomer density is different from the actual density. Moreover, if the coefficient $\gamma$ of $\gamma S_{1}$ in Eq. (12) is taken to be slightly larger than 1, the volume is expected to be slightly smaller than the current $V$; consequently, $f$ becomes smaller and is hence controllable.

The snapshots of the PVDF model are shown in Figs. 2(a) and (b), where the external electric field is $E\!=\!0$ and $E\!=\!0.15$, corresponding to $E_{\rm exp}\!\simeq\!189$ MV/m. The small red cylinders on the surface denote the variable $\sigma$. $\sigma$ aligns along the $z$ axis even when $E\!=\!0$ because $\lambda$ is set as $\lambda\!=\!0.345$, corresponding to the ferroelectric behavior of the $\beta$ phase of the PVDF.

According to Fig. 3(a), the simulation results of $P$ vs. $E$ are in agreement with $P_{\rm exp}$ vs. $E_{\rm exp}$. The influence of the other simulation parameters $\kappa$ and $v_{0}$, which are set as $\kappa\!=\!1.8$ and $v_{0}\!=\!0.1$, on the PE curve is almost negligible.

Fig. 3(b) shows the relation between the height strain $\varepsilon_{H}$ and electric field. 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\ell_{0}$, where $\ell_{0}$ is the mean bond length. We should note that the experimental loop has a 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 decreases 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 exhibits almost no hysteresis [4]. This finding implies that the strain is directly related to the polarization rather than the electric field. Therefore, we expect that the electrostriction constant, i.e., the slope of the strain-polarization curve, can be controlled by $\kappa$. The solid line in Fig. 4(a), denoted “Theory”, is plotted using $P^{2}({S})\!=\!P^{2}(0)\!+\!(1/\kappa_{33}){S}$ with the electrostriction constant $\kappa_{33}\!=\!-2.4$ and experimental remnant polarization $P(0)$, where $S$ denotes the strain. The constant $|\kappa_{33}|$ decreases with increasing $\kappa$. The experimental sample involves the constant $\kappa_{33}\!=\!-2.4$, and this value is achieved in our model as $\kappa\!=\!1.8$.

The elongation of the cylinder diameter $D$ is calculated as $D\!=\!2\sqrt{V/(\pi H_{0})}$ with the initial height $H_{0}$ at $E\!=\!0$. Using this $D$, we obtain the diameter strain $\varepsilon_{D}$, which is considered to be the “nominal strain” because $H_{0}$ is used to calculate $D$; hence, we estimate the Poisson ratio $\nu\!=\!\!-\varepsilon_{D}/\varepsilon_{H}$ (Fig. 4(b)). $\nu$ tends to be close to 0.3 with increasing strain within the experimental electric field range. This value is comparable to the experimental value of 0.33 for $\beta$-PVDF obtained through mechanical testing [24].

Thus, we conclude that the FG model successfully reproduces the electromechanical properties of $\beta$-PVDF. Notably, the PE, strain-electric (SE) and Poisson ratio data are obtained with a single set of parameters.

To ensure that the simulations are correctly performed, we plot $S_{1}/N$ vs. $E$, as shown in Fig. 5(a), where the simulation unit is used for all the quantities including $E$. As mentioned above, $S_{1}/N\!=\!1.5$ is satisfied. This finding implies that the constraints imposed on the lattice, such as the hard wall, do not influence the equilibrium property because this relation is satisfied only when no constraint is imposed on the vertex positions ${\vec{r}}_{i}$ [13, 14, 15]. The relation of volume $V$ and $E$ in Fig. 5(b) changes symmetrically with respect to the changes in $E$ and $-E$, as expected.

We describe (i) the advantage of the FG model over the other models and (ii) how the FG model facilitates the understanding of such unusual and complex phenomena. To highlight these points, we intuitively express the basic idea. First, we must emphasize that the elastic energy or bond potential $\gamma S_{1}\!=\!\gamma\sum_{ij}\Gamma_{ij}l_{ij}^{2}$ in Eq. (8) is of the form of the sum of the “effective coupling constant” $\gamma\Gamma_{ij}$ and squared length $\ell_{ij}^{2}$ between molecules. The term “molecules” does not always express real molecules and is used in the context of coarse-grained particles, as described above.

In the original model, in which the Euclidean metric is assumed, the effective coupling constant is $\gamma\Gamma_{ij}\!=\!\gamma(\Leftrightarrow\Gamma_{ij}\!=\!1)$, and therefore, the elongation of the materials governed by this bond potential must be isotropic. However, the elongation of the polymers under uniaxial forces, such as that of PVDF under electric fields, is anisotropic or direction dependent, and therefore, the former part, $\lambda$, of the bond potential should be direction dependent, because the latter part $\ell_{ij}^{2}$ is position dependent and does not depend on the polarization direction.

The question remains regarding the dependency of the constant $\lambda$ on the direction. One potential solution is to consider that the metric for the length unit of the local coordinates inside the material is dependent on the direction of the polarization vectors. In this manner, the deformation, such as that of PVDF, can be understood. The bond potential of the FG model is dependent on not only the squares of the distance but also the coupling constant, as described above. Moreover, the distance increases (decreases) for a small (large) coupling constant, as mentioned in Section 2, because the mean value of $S_{1}$ remains unchanged from that of the original model even if $S_{1}$ is deformed by the FG modeling prescription. Accordingly, the coupling constant ($\Leftrightarrow\Gamma_{ij}$) can be numerically or automatically determined by performing averaging over the scales of the molecular distance altered by fluctuations in the polarization direction.

The key aspect is that the macroscopic property, such as the elongation of PVDF, varies depending on a certain position-dependent microscopic physical quantity of the material, such as the polarization vector of the PVDF. Depending on the value of the microscopic quantity, the unit Finsler length $v_{ij}$ is suitably defined, similar to that in Eq. (2) with a suitable value of $v_{0}$, and the effective coupling constant ($\Leftrightarrow\Gamma_{ij}$) is automatically determined. In general, the macroscopic property obtained in this manner in the FG modeling is independent of the directional degree of freedom of the microscopic quantity and detailed information of the physical process such as the interactions of the electrons and atoms [25]. Consequently, the FG modeling technique can be implemented easily in many models for complex phenomena and can be applied in other modeling techniques.

We demonstrate that the FG model can be applied to examine the electromechanical properties of ferroelectric polymers. The Monte Carlo simulation results are in agreement with the experimental data of $\beta$-PVDF. Both PE and SE field curves are reproduced using a single set of simulation parameters. Moreover, the simulated Poisson ratio is reasonable. The technique used for the PVDF is also applicable to ferroelectric ceramics such as ${\rm BaTiO_{3}}$. The magnetostriction of ferromagnetic materials will be examined in future work.

The author H.K. thanks Laurent Chazeau for the valuable discussions on hysteresis curves. This work is supported in part by the Collaborative Research Project of the Institute of Fluid Science, Tohoku University, and by a JSPS Grant-in-Aid for Scientific Research on Innovative Areas ”Discrete Geometric Analysis for Materials Design”: Grant Number 20H04647 and JSPS KAKENHI Grant Number JP17K05149. The author V.E. thanks President Dr. Hiroshi Fukumura of Sendai KOSEN for the warm hospitality during the four-month stay from 2019 to 2020, which was supported in part by JSPS KAKENHI Grant Number JP17K05149.