A combination of experimental and numerical method for optimizing the sensitivity of ultra-thin piezoelectric sensor with interdigitated electrodes

We report the development of a finite element model (FE-model) for ultra-thin piezoelectric poly(vinylidene-trifluoroethylene) sensor with interdigitated electrodes (IDE) which includes the effect of a non-homogenous poling field determined via the combination of experimental and numerical methods. The non-homogenous poling magnitude is estimated by comparing the remanent polarization (P r) of IDE based device to the P r of the same material in metal–insulator–metal electrode configuration. The non-homogenous poling orientation is estimated by comparing the experimentally determined normal mode sensitivity (S n) values to FE-modelled sensitivity values with different poling orientation distributions. The poling orientation distribution is modelled using two approaches: (a) 33-direction parallel and perpendicular to the electrode plane and (b) 33-direction defined by an average angle. The first approach yields the best correspondence with the experimental results (R 2 = 94.70% and σ = 0.10 pC N−1) and it is used to optimize the device geometry and poling condition for maximum S n.


Introduction
Piezoelectric polymeric materials based on poly(vinylidenefluoride-trifluoroethylene) P(VDF-TrFE) have been recently employed in the fabrication of ultra-thin sensors which enable user-friendly and comfortable measurement of mechanically induced bio-signals [1][2][3][4]. However, traditional piezoelectric sensors with metal-insulator-metal (MIM) electrode configuration have limitations in this regard because decreasing their thickness increases their capacitance and thereby results in lower voltage sensitivity of the device. Novel electrode configurations based on interdigitated electrode (IDE) structure have been proposed to combat this issue; in this case, the device capacitance does not depend on the thickness of the piezoelectric layer to the same extent (see e.g. [5] and table S1), and therefore the device voltage sensitivity is less affected by the decreasing piezoelectric layer thickness. This enables the fabrication of very thin, yet high sensitivity devices as demonstrated by Lozano et al [1].
However, in addition to the capacitance, the device sensitivity also depends on the charge generation ability of the piezoelectric material, which for MIM based devices is well understood, but for IDE based devices less so: the main challenge in understanding the IDE based device charge generation is related to the non-uniform poling field between the consecutive electrodes. As the direction of the electric field lines determines the poling direction of the material, for example the positive 33-direction in IDE based device may point parallel or perpendicular to the electrode plane, or any angle between them, depending on the point under investigation (i.e. various 33-directions will encompass all of 360 • ; see figure 1(A)); compare this to a MIM based device where the 33-direction always points to one direction i.e. perpendicular to the electrode plane. Furthermore, the magnitude of the poling varies significantly in the case of IDE, whereas for MIM it is homogenous throughout the material: electrostatic FE-modelling [6][7][8][9][10] has shown that the IDE-based devices have an inactive area close to the centerline of the electrode where the electric field magnitude, and therefore also the poling magnitude, is zero. Thus, to fully understand the factors affecting the sensitivity of the IDE based device, one would have to solve for the net effect generated by the various magnitudes and directions of poling throughout the material. Beckert and Kreher [6] used a two-step approach to achieve this: they first modelled the non-uniform electric field and then used these results to assign a local poling direction/magnitude for each element of the piezoelectric model. However, this may be challenging because of the time and computational power required in dividing the piezoelectric material model into small elements, assigning each of these with a local poling direction/magnitude and finally solving the complex model. Toprak and Tigli [9] used similar two step approach; however, they simplified the latter step by using an average poling orientation/magnitude on top of the electrode instead of element specific local orientation/magnitude.
In this article we try to improve upon the latter methods by basing the determination of nonhomogenous poling orientation/magnitude on experimental results instead of electrostatic FEmodelling. The non-homogenous poling magnitude is estimated based on remanent polarization (P r ) measurement. Specifically, we assume that the poling magnitude in the active area is 100%, and in the inactive area 0%. The width of the active and inactive area is then estimated by comparing the measured P r value for the IDE to the P r of the same material in the MIM configuration. The non-homogenous poling orientation is estimated by comparing the experimentally determined sensitivity values to FEmodelled sensitivity values with different poling orientation distributions. Specifically, the poling orientation distribution is included in the FE-model using two approaches: (a) the active area is divided into poling orientations which are parallel and perpendicular to the electrode plane (see figure 1(B)) and (b) the poling orientation in the active area is described by average angle as suggested by Toprak and Tigli [9] (see figure 1(C)). A statistical analysis is then performed to predict the poling orientation distribution based on the sample dimensions and the poling condition, and the resulting regression equation is used as an input in the FE-model. We show that the approach (a) yields the best match between the FEmodelled and experimentally determined sensitivity values. Finally, we use the thus generated FE-model to optimize the sample dimensions and the poling condition for maximum sensitivity.

Results and discussion
2.1. Sample description IDE based piezoelectric samples with inkjet printed PEDOT:PSS interdigitated electrodes and bar coated P(VDF-TrFE) layer were fabricated on top of a Parylene-C polymer as described in our previous publication [1]. Seventeen IDE samples were included in the analysis with P(VDF-TrFE) thickness (t P(VDF-TrFE) ) ranging from 5.7 µm to 25.2 µm and electrode width (w) ranging from 134 µm to 156 µm; the electrode-to-electrode gap (g) is dependent on electrode pitch (200 µm) and ranges from 44 µm to 66 µm (see figure 1). The normal mode sensitivity (S n ) and remanent polarization (P r ) of each IDE sample was measured as described in experimental section. The results are summarized in table S1.

Non-homogenous poling magnitude
Remanent polarization (P r ) values extracted from the polarization-electric field loop (PE-loop) measurements were employed in determining the width of the inactive area at the center of the finger electrodes. The remanent polarization of the IDE samples (P r.IDE ) is extracted from PE-loop measurement where the polarization (P) is determined as (see figure 2): where Q is the measured charge, A the IDE electrode area from which the charge is measured, N the number of electrodes, L the overlapping length of the electrodes and w the width of the electrodes [1,11,12]. However, this does not take into account the inactive area in the middle of the finger electrodes. A more realistic approach would be to assume that the P(VDF-TrFE) reaches its maximum poling magnitude (P r = 7.3 µC cm −2 as measured using the MIM reference device, see figure 2) in the active area while the poling magnitude is assumed to be 0 µC cm −2 in the inactive area. In this case the equation (1) can be written as where P active is the polarization measured from the active area, A active the active area and w active the width of the active area. Because the measured charge must be the same in both cases, we end up with the following equation: As per our assumption, we may now replace the P active with MIM remanent polarization (P r.MIM = 7.3 µC cm −2 , see figure 2) while the P is the P r.IDE extracted from the original PE-loop measurement: The width of the inactive area (w inactive ) is then simply: The results are summarized in table S1.

Non-homogenous poling orientation
A 3D FE-model for the IDE based piezoelectric device in the S n measurement situation was generated using COMSOL software with piezoelectric multiphysics module. The geometry was miniaturized such that it consisted of the representative volume element (RVE) of the IDE (see figures 1(B)-(D)). The inactive area was omitted from the model as it is assumed to have 0% poling magnitude and it does not therefore contribute toward piezoelectric charge generation. In contrast, the poling magnitude in the active area is assumed to be complete (100%) and the piezoelectric properties are therefore described by the full coupling matrix of P(VDF-TrFE) [13]. These and other relevant material parameters are given in SI. The effect of non-homogenous poling orientation was included in the model using two approaches: Model 1) The active area on top of the electrode is divided to areas where the 33-direction is perpendicular (in negative/positive z-direction, see figure 1(B)) and parallel to the electrode plane (in the positive y-direction). A unitless parameter z poling , with values between 0 and 1, was used to describe the poling direction distribution; for example, z poling = 0.4 indicates that 40% of the active area on top the electrode is poled in direction perpendicular to electrode plane and 60% in direction parallel to electrode plane. Model 2) The poling direction is uniform over the active area of the electrode, but the 33-direction is at an angle when compared to the z-axis (see figure 1(C)). This corresponds to the average poling direction approach suggested by Toprak and Tigli [9]. z angle is used to describe the angle between the 33-direction and z-axis and it can take values between 0 and 90 • .
Finally, a compressive force was applied in negative z-direction with pressure equal to the experimental situation while evaluating the surface charge generated on one of the electrodes. This was translated to total charge output (Q out ) by multiplying it with the number of RVEs activated in the S n measurement: where N RVE is the number of RVEs, Q s the surface charge, D n the normal component of the electric displacement field and S the surface where the D n is evaluated (i.e. the electrode surface). Considering the RVE geometry where the width and depth of the electrode is equal to 0.5·w active , we can write: However, the N RVE decreases as the w active increases. Consider a sample where the piezoelectric material on top of N electrode pairs with overlapping length L is compressed uniformly. Then: Combining equations (4)-(8) yields: The S n can be then solved using the formula: where F is the applied compressive force. It is important to note that the relation S n ∝ P r.ratio ·w is valid for all force probe geometries (not just the rectangular probe assumed above).
The following stepwise approach was used to relate the poling orientation distribution to the sample geometry and poling condition: (a) FE-models for each sample were solved for S n while sweeping the 1. z poling from 0.1 to 0.9 2. z angle from 0 to 90 • (b) Regression lines were fitted to the modelled S n data and solved for the z poling and z angle values which produce the measured S n values. (c) The P r.ratio , t P(VDF-TrFE) and w were used to create regression model for the z poling and z angle .
The results of steps 1 and 2 are plotted in figures 3(A) and (B). For all thicknesses, the sensitivity increases as the active area with the poling direction perpendicular to the electrode plane increases (i.e. increase in z poling or decrease in z angle ). This can be qualitatively explained by considering the direction of dominating forces and the magnitude of the activated piezoelectric coefficients in the different poling direction areas of the material: the dominating forces in the S n measurement are perpendicular to the electrode plane and this leads to activation of d 33 -coefficient (−34.9 pC N −1 ) in the perpendicularly poled material, and activation of lower value d-coefficients (e.g. d 31 -coefficient 10.4 pC N −1 ) in the parallel poled material. Thus, increasing the portion of perpendicularly poled material should lead to higher sensitivity. The z poling and z angle values producing the experimentally determined S n value can be then solved from the regression equation as shown in figures 3(C) and (D). The solved z poling and z angle values are summarized in table 1.
The sample geometry and poling condition can be then related to the solved poling orientation by creating a regression model based on the analysis of variance (ANOVA) method. The predictors included in the ANOVA were the P r.ratio , t P(VDF-TrFE) and w. Table S3 shows the ANOVA results for the z poling after removing the insignificant predictors using backwards elimination with P-value = 0.05 as a criterion; and z angle = 653 − 10.833t P(VDF−TrFE) − 439.6P r.ratio + 11.16w + 0.08754t P(VDF−TrFE) · t P(VDF−TrFE) + 330.3P r.ratio · P r.ratio − 0.0399w · w where t P(VDF-TrFE) and w are expressed in µm. Both models fit the data well as is indicated by the high coefficient of determination (R 2 ) value and low standard deviation of residual (σ): for z poling R 2 = 98.69% and σ = 0.0201, and for z angle R 2 = 99.52% and σ = 0.9265. Based on the F-value indicated in tables S3 and S4, the thickness has the strongest effect followed by quadratic term on the P r.ratio , P r.ratio , and w. Graphical interpretation of these results is shown in figure 4 for z poling : increasing the thickness of the P(VDF-TrFE) tends to decrease the amount of poling perpendicular to the electrode plane; P r.ratio has a similar effect with minimum z poling ≈ 0.75, after which the perpendicular poling starts to again increase; the electrode width has the opposite trend with maximum z poling around 148 µm. Figure 5 shows that the larger electrode widths decrease the effect of thickness on the z poling . A similar non-linear relationship between the piezoelectric layer thickness, the electrode width and the poling orientation has been previously observed by e.g. Toprak and Tigli [9].
The regression equations for z poling (equation (11)) and z angle (equation (12)) can be then used as an  .56% for z poling and z angle based FEmodels, respectively. The z poling based FE-model is therefore the preferable approach for optimization of the geometry and poling condition of the IDE based piezoelectric sensor.

Optimizing the geometry and poling condition for maximum S n
Now that the correspondence between the z poling based FEM and experiments has been verified, the FEM can be used to optimize the geometry of the device to reach the maximum S n value. A design of experiments (DoE) approach was adopted to this end with P r.ratio , t P(VDF-TrFE) and w as input and FEmodelled S n as output; a full factorial central composite design with 20 different parameter combinations  was chosen as this allows the estimation of possible quadratic effects. These are expected due to the previously observed quadratic dependence of z poling on P r.ratio , t P(VDF-TrFE) and w. The DoE runs and resulting FEM output are shown in tables S5 and S6 shows the ANOVA results after backwards elimination of insignificant predictors using P-value = 0.05 as criterion.
The regression model generated based on the ANOVA results is: − 45.89P r.ratio − 0.00792w · w + 34.90P r.ratio · P r.ratio where t P(VDF-TrFE) and w are in units µm and S n is in units pC N −1 . The model fits the data well as is indicated by the high R 2 value of 96.85% and low σ value of 0.1863. The F-values shown in the table S6 indicate that the quadratic term of P r.ratio has the strongest effect followed by t P(VDF-TrFE) and w. A graphical interpretation of the results is shown in figure 6. The normal mode sensitivity decreases linearly as the thickness of P(VDF-TrFE) is increased, and the P r.ratio and w have minimum and maximum values occurring approximately in the middle of the design space. The parabolic shape originates from the non-linear dependence of poling direction distribution on P r.ratio and w. Omitting the effect of poling direction distribution, increasing P r.ratio would be expected to increase the S n linearly because as the P r.ratio is increased, the piezoelectrically active width (see equation (10)) increases linearly. The same applies for the electrode width. However, in both cases the non-linear effect induced by the poling direction distribution complicates the shape of the observed effects. The interactive effect of t P(VDF-TrFE) and w (see figure 7) shows that the effect of t P(VDF-TrFE) on S n decreases as electrode width is increased.
It is interesting to note that the normal mode sensitivity can be increased by reducing the P(VDF-TrFE) thickness. However, we also observed a significant drop in the yield of the poling step due to electric arcing between consecutive IDEs once the P(VDF-TrFE) thickness was reduced below ∼5 µm. This practical limitation was taken into account when determining the boundary conditions for the following parameter optimization.
Finally, the equation (13) can be used to find the optimum sample geometry and poling condition which results in maximum S n value. The ten best solutions are ranked according to the composite desirability in table S7; the optimum parameter combination occurs at minimum t P(VDF-TrFE) (5 µm), intermediate w (∼144 µm) and maximum P r.ratio (0.9).

Conclusion
A combination of experimental results, statistical modelling and FE-modelling was used to model the non-homogenous poling field of an ultra-thin IDE based piezoelectric sensor. The non-homogenous poling magnitude was estimated based on remanent polarization measurements and two FE-modelling approaches for the non-homogenous poling orientation were compared. The FE-model with the best approach had good correspondence with the experimental results (R 2 = 94.70% and σ = 0.10 pC N −1 ) and it was used to optimize the device geometry and poling condition for maximum normal mode sensitivity.

Fabrication of the P(VDF-TrFE) sensor
The IDE based and MIM based reference samples were fabricated as described in our previous publications [1] and [2], respectively. See figure S1 for representative microscope images of the IDEs, piezoelectric layer and a final sample.

Sensor dimensions
The thickness of the P(VDF-TrFE) layer was measured using a stylus profilometer (Dektak XT, Bruker). Each sample was scanned from two locations and the average of these two measurements was used as the thickness. The IDE finger width was measured with an optical microscope (BX60M, Olympus) from three locations and the average of these measurements was calculated.

Ferroelectric, piezoelectric and electric characterization
The PE loops were measured using the TF 2000 Analyzer (AIXACCT) coupled to a high voltage amplifier (610C, TREK). For IDE based devices, the applied electric field (E) was calculated as the ratio of the applied voltage (V) and the gap (g) between the consecutive electrodes i.e. E = V/G; for MIM based device it was calculated by as ratio of voltage and the thickness of the P(VDF-TrFE) layer (t) i.e. E = V/t. The remanent polarization (P r ) values were extracted from the saturated PE-loops.
The normal mode sensitivity (S n ) was measured using PiezoMeter PM300 (Piezotest Ltd) with flat circular probe with a radius of 5 mm while applying 0.25 N of dynamic force with 110 Hz frequency; the sample was clamped between the probes using 10 N of static force.
The capacitance values were measured at 1 kHz frequency using a semiconductor analyzer (B1500A, Keysight).

Modelling
COMSOL Multiphysics (v.5.6.) with piezoelectric module was used to generate the FE-models. Origin (v. 2019b) was used to fit the regression lines to the FE-modelled data. Minitab Statistical Software (v. 21) was used to perform ANOVA for z poling and z angle results; generate and analyze the DoE; and find the optimum parameters for maximum S n .

Data availability statement
The data generated and/or analyzed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.