Domain splitting and enhanced piezoelectricity in ferroelectric thin films with square grid structure

The performances of ferroelectric (FE) thin films and nano-structures have been receiving considerable interest not only due to the significance in fundamental researches but also for their diverse applications in piezoelectric micro-systems, FE memories, high-frequency electronic components, etc. [1–3]. Those applications rely substantially on the FE polarization and electromechanical (piezoelectric) behaviors, which are related to the FE domain structures and their responses to different environmental stimuli such as electrical field, stress, spatial confinement (dimensionality), etc. Along this line, extensive investigations on the FE domain structure and its control in order to optimize these properties have been reported [4–6]. The so-called domain engineering on ferroelectrics has become of high concern in recent years since those application performances can be then remarkably improved [7,8].

The most representative class of FE materials targeting such applications includes FE perovskite oxides like $\text{BaTiO}_{3}$ and $\text{PbTiO}_{3}$ which have tetragonal FE phase originating from cubic paraelectric phase upon decreasing temperature T [9,10]. For a FE system like $\text{BaTiO}_{3}$, the well-aligned stripe-like $90^{\circ}$ domain structure is one of the core characteristics [9]. For the sake of simplification, we discuss a two-dimensional (2D) FE lattice as an approach to FE thin film, where the FE polarizations are confined onto the xy-plane along the [01] and [10] easy axes. This $90^{\circ}$ domain structure is the consequence of the delicate balance mainly between the electric dipole-dipole interactions (depolarization effect) associated with the FE polarization and the electromechanical (elastic plus electrostrictive) energy associated with the piezoelectricity and electrostriction effect, while the domain wall energy may not be negligible, too [11]. This delicate balance is the reason for sensitive responses of the domain structure to either external electric field or stress/pressure, so as to exhibit excellent electromechanical performances in addition to the favored FE properties [12].

It is noted that both the dipole-dipole interactions and the electromechanical energy in a FE lattice are long-ranged and they can be strongly modulated by imposing various spatial constraints, in particular the dimensionality and geometric constraints. They are responsible for very scattered FE and piezoelectric properties of finite-size FE materials such as nano- or micro-scale dots and islands as well as substrates confined thin films [13–15]. Consequently, the regular stripe-like or twinned $90^{\circ}$ domain structure in bulk lattice may be replaced by vortex-like domain pattern and others. Quite a number of works reported observations of various fascinating domain patterns by modulating these spatial constraints [16,17]. Taking a naked FE island (in-plane size L and thickness D, and $D/L\ll 1$ for thin films) with free boundary as an example, one can easily predict the evolution of domain structure, with decreasing L from stripe-like domain pattern into a coexistence of the stripe-like pattern and the vertex-like one, eventually into a single vertex pattern [18]. For such small-size island, the dipole-dipole interaction in competing with the electromechanical energy dominates the formation of vertex structure. It was reported that such FE islands in some cases do exhibit better electromechanical performances than their bulk counterpart.

Given the rapid progress in micro-fabrication technology of materials, one may have additional approaches to engineer the FE domain structure. While a fabrication of well-aligned FE island lattices (arrays) faces technical challenging, it represents one end of releasing the spatial constraints. In this case, the dipole-dipole interaction and electromechanical energy are seriously released. The other end would be an alternative scheme in which an infinite FE thin film can be patterned by removing finite-size units, so that a grid thin film structure with periodically aligned voids (void array) can be fabricated. This end scheme seems quite less touched in the literature so far but offers several particulars. Starting from a stripe-like $90^{\circ}$ domain structure, fabricating such voids would result in several consequences. First, the thin film as a whole remains continuous, but the voids will partially release the electromechanical energy of the neighboring film region while the dipole-dipole interaction enforces the electric polarizations surrounding the voids to re-align. By properly designing the spatial distribution, shape, and size of the voids, one is allowed to modulate substantially the stripe-like domain pattern for specific purposes. Second, due to the high mobility of the $90^{\circ}$ domain walls, such voids may benefit to the enhancement of the FE and electromechanical responses. Third, with respect to the FE island array, this grid thin film structure can be easily and efficiently fabricated by up-bottom fabrication techniques, which is technologically favored and cost-competitive.

In this work, we investigate the spatial evolution of the $90^{\circ}$ domain structure in a FE thin film with periodically aligned square void array, based on the well-known Landau-Devonshire phenomenological theory [19–22]. Certainly, the shape of voids can be other than square but the underlying physics should be similar. We start from a 2D $L \times L$ square lattice with periodic boundary conditions as an approach to an infinite FE thin film. On each site i two order parameters are imposed: polarization vector $\boldsymbol{P}(\boldsymbol{r})=(P_{x}$, $P_{y})$ and displacement vector $\boldsymbol{u}(\boldsymbol{r})= (u_{x}$, $u_{y})$, where x and y are the coordinate axes. Following earlier literature, the free energy of the FE lattice includes the Landau energy, gradient energy, dipole-dipole interaction, elastic energy, electrostrictive interaction, and electrostatic energy [18]. For such a lattice, the total free energy density can be written as

$$\begin{equation} f_{\textit{total}} =f_{ld} +f_{g} +f_{\textit{dip}} +f_{\textit{ela}} +f_{es} +f_{E}, \end{equation} \tag{ 1 }$$

where the terms $f_{ld}$, $f_{g}$, $f_{\textit{dip}}$, $f_{\textit{ela}}$, $f_{es}$, $f_{E}$ are the Landau energy, gradient energy, dipole-dipole interaction, elastic energy, electrostrictive energy, and electrostatic energy, respectively. The electromechanical energy includes $f_{\textit{ela}}$ plus $f_{es}$. The Landau energy is a polynomial extended to the sixth order:

$$\begin{eqnarray} f_{ld} &=A_{1} \left({P_{x}^{2} +P_{y}^{2} } \right)+A_{11} \left({P_{x}^{4} +P_{y}^{4} } \right)+A_{12} \left({P_{x}^{2} P_{y}^{2} } \right)\nonumber\\ &\quad +\,A_{111} \left( {P_{x}^{6} +P_{y}^{6} } \right)+A_{112} \left( {P_{x}^{4} P_{y}^{2} +P_{y}^{4} P_{x}^{2}}\right), \end{eqnarray} \tag{ 2 }$$

where A with various subscripts are the Landau coefficients. Usually, $A_{1}$ is T-dependent and can be written as $A_{1}=A_{10}(T-T_{0})$, where $T_{0}$ is the critical point for the FE transition, and $A_{10}= (2\varepsilon_{0}C)^{-1}$ with vacuum permittivity $\varepsilon_{0}$ and Curie constant C. Positive $A_{1}$ and negative $A_{11}$ are required for the FE transitions. The domain wall energy is expressed by the gradient energy form:

$$\begin{eqnarray} f_{g} &=\frac{1}{2}\Biggl[G_{11} \left({P_{x,x}^{2} +P_{y,y}^{2}}\right)+G_{12} P_{x,x} P_{y,y}\nonumber\\ &\quad +\,G_{44} \left( {P_{x,y} -P_{y,x}}\right)^{2}+G_{44}^{\prime} \left( {P_{x,y} +P_{y,x} } \right)^{2}\Biggr], \end{eqnarray} \tag{ 3 }$$

where $P_{i,j}=\partial P_{i}/\partial r_{j}$, and $G_{ij}>0$ are the gradient energy coefficients. The long-range electric dipole-dipole interaction is written as

$$\begin{eqnarray} f_{\textit{dip}} &=\frac{1}{4\pi \varepsilon }\sum\limits_{\left\langle j \right\rangle}\Biggl[\frac{P_{i} \cdot P_{j} }{\left| {r_{i} -r_{j} } \right|^{3}}\nonumber\\ &\quad -\,\frac{3P_{i} \cdot \left( {r_{i} -r_{j} } \right)\left[ {P_{j} \cdot \left( {r_{i} -r_{j} } \right)} \right]}{\left| {r_{i} -r_{j} } \right|^{5}}\Biggr],\quad \end{eqnarray} \tag{ 4 }$$

where the pre-factor $(1/4\pi \varepsilon )=1$ is taken for normalizing purpose. We can re-write eq. (4) as

$$\begin{equation} f_{\textit{dip}} =-\frac{1}{2}E_{1} \cdot P_{i}, \end{equation} \tag{ 5 }$$

where the internal electric field $E_{1}$ is

$$\begin{equation} E_{1} =-\frac{1}{2\pi \varepsilon }\sum\limits_{\left\langle j \right\rangle } {\left[ {\frac{P_{j} }{\left| {r_{i} -r_{j} } \right|^{3}}-\frac{3\left( {r_{i} -r_{j} } \right)\left[ {P_{j} \cdot \left( {r_{i} -r_{j} } \right)} \right]}{\left| {r_{i} -r_{j} } \right|^{5}}} \right]}. \end{equation} \tag{ 6 }$$

Due to the existence of voids in the lattice, big lattice size L has to be taken. A precise account of this long-range term using the Ewald summation scheme becomes practically overtedious [23]. Fortunately, for 2D lattice, one can choose a preset truncating distance $R_{0}$ from the central site in order make the computation tractable. This scheme is sufficient for a reliable account of this interaction ($R_{0}=8$ in our simulation) [24,25].

Subsequently, the elastic energy reads

$$\begin{equation} f_{\textit{ela}} =\frac{1}{2}C_{11} \left( {u_{x,x}^{2} +u_{y,y}^{2} } \right)+C_{12} u_{x,x} u_{y,y} +\frac{1}{2}C_{44} u_{x,y}^{2}, \end{equation} \tag{ 7 }$$

where $u_{i,i}=\partial u_{i}/\partial r_{i}$, $u_{i,j}=\partial u_{i}/\partial r_{j}+\partial u_{j}/\partial r_{i}$, and $C_{ij}$ are the elastic coefficients. The electrostrictive energy accounting the coupling between the polarization field and strain field is written as

$$\begin{equation} f_{es} =-u_{x,x} \eta_{xx} -u_{y,y} \eta_{yy} -u_{x,y} \eta_{xy}, \end{equation} \tag{ 8 }$$

with strain $\eta_{xx}=Q_{11}P_{x}^{2}+Q_{12}P_{y}^{2}$, $\eta _{yy}=Q_{11}P_{y}^{2}+Q_{12}P_{x}^{2}$ and $\eta_{xy}=Q_{44}P_{x}P_{y}$, where $Q_{44}$ are the electrostrictive coefficients. Finally, the electrostatic energy is

$$\begin{equation} f_{E} =-E\cdot P_{i}, \end{equation} \tag{ 9 }$$

with applied electric field $E = (E_{x}, E_{y})$.

The void array is simply introduced into the FE lattice by imposing $\boldsymbol{P}(\boldsymbol{r}) \equiv 0$ and $u_{i,j}(\boldsymbol{r})\equiv 0\ (i, j=x, y)$ in the pre-designed regions. Here we only consider square voids with the two sides aligning along the x-axis and y-axis, respectively. It is noted that these voids in the lattice makes it mathematically challenging to properly handle the mechanical balance conditions on the boundaries between thin film matrix and voids. This difficulty becomes particularly clear if the time-dependent Ginzburg-Landau equation (TDGL, Langevin-type) is used to solve this issue numerically [16]. However, no such difficulty appears if one employs the Monte Carlo (MC) simulation to track the evolution of the FE lattice [18]. A reasonable assumption here is that the response time of elastic strain is far shorter than that for dipole relaxation, such that the mechanical equilibrium is instantly established when an electric dipole flip event occurs [15,26]. The object lattice of our simulation is $\text{BaTiO}_{3}$ thin film. The relevant materials parameters are taken from the literature [18]. We assume that the lattice constant $a_{0}$, Landau coefficient $\alpha_{0}$, and polarization $P_{0}$ at $T=T_{1}$ far below $T_{0}$ are known as reference, where $\alpha_{0}=\alpha_{10}(T_{1}-T_{0})$ and $\alpha_{10} >0$. All the free energy terms are normalized by term $\vert\alpha_{0}\vert P_{0}^{2}$ which has an energy unit [27]. Table 1 lists the chosen parameters for simulation.

Table 1:.  Parameters chosen for the simulation ($1/4\pi \varepsilon =1$, Boltzmann constant $k_{B}=1$).

Parameter (unit)	Value	Parameter (unit)	Value	Parameter (unit)	Value
$A_{10}\ (\alpha_{10})$	0.10	$A_{11}\ (\alpha_{11}P_0^2/ \vert\alpha_{0}\vert)$	−0.24	$A_{12}\ (\alpha_{12}P_0^2/\vert\alpha_{0}\vert)$	4.50
$A_{111}\ (\alpha_{111}P_0^4/\vert\alpha_{0}\vert)$	0.49	$A_{112}\ (\alpha_{112}P_0^4 /\vert \alpha_{0}\vert)$	1.20	$G_{11}\ (g_{11}/a_0^2 \vert\alpha_{0}\vert)$	1.60
$G_{12}\ (g_{12}/a_0^2 \vert \alpha_{0}\vert)$	0.00	$G_{44}\ (g_{44}/a_0^2 \vert\alpha_{0}\vert)$	0.80	$G^{\prime}_{44}\ (g^{\prime}_{44}/a_0^2 \vert\alpha_{0}\vert)$	0.80
$C_{11}\ (c_{11}/\vert\alpha_{0}\vert P_0^2)$	2.75	$C_{12}\ (c_{12}/\vert\alpha_{0}\vert P_0^2)$	1.79	$C_{44}\ (c_{44}/\vert\alpha_{0}\vert P_0^2)$	0.54
$Q_{11}\ (q_{11}/\vert\alpha_{0}\vert)$	0.142	$Q_{12}\ (q_{12}/\vert\alpha_{0}\vert)$	−0.0074	$Q_{44}\ (q_{44}/\vert\alpha_{0}\vert)$	0.0157
$L(a_{0})$	64	$T_{0}(T_{c}/(T_{c}\text{-}T_{1}))$	5.00

For practical simulation, we submit a FE lattice free of any void to a high $T \gg T_{0}$ for annealing. Then the lattice is gradually cooled with sufficient relaxation time down to a low $T \ll T_{0}$. The details of the simulation were reported in our earlier work [24,25] and no description is given here. Figure 1(a) shows the simulated $90^{\circ}$ domain structure with L = 64. The well-aligned stripe-like $90^{\circ}$ domain walls along the $[1\bar{1}]$ direction are observed. The results are well consistent with phase-field simulations using the TDGL equation [15,26,28]. Subsequently, a periodically distributed square void array of given size $L_{s}$ is introduced into the lattice and a sufficiently long-time relaxation of the domain structure is allowed. The eventually reached domain structure as the energy lowest state (the ground state) is checked carefully by choosing various annealing paths.

Zoom In Zoom Out Reset image size

We perform extensive simulations via the standard Metropolis algorithm on the domain structure evolution by imposing different spatial distribution and size of the voids. It is shown that the evolutions share qualitatively similar characteristics. As an example, we present in figs. 1(b)–(d) the simulated domain structures with a $4\times 4$ void array, forming a grid FE thin film. The void size $L_{s}$ takes 0, 4, 6, and 12, respectively. At the first glance, the existence of voids in the FE thin film has marked impact on the domain structure. However, differently from FE islands where the vertex pattern is preferred, we observe an interesting pattern evolution when $L_{s}$ takes different values. At $L_{s}\le 4$, the well-aligned $90^{\circ}$ domain structure with straight walls (fig. 1(a)) can be maintained roughly, whereas only local deformation in the domain walls neighboring the void boundaries is seen, as shown in fig. 1(b). This local deformation becomes complicated as $L_{s}>4$ and the pattern at $L_{s}=6$ is shown in fig. 1(c). Surprisingly, regular domain splitting sequence is clearly observed due to the introduction of the void array. By such splitting, one can halve the domain width and periodicity. This domain splitting can sustain until $L_{s}$ is big enough, such as $L_{s}=12$, as shown in fig. 1(d), where the FE lattice consists of crossing thin grids along the x-axis and y-axis. Due to the strong depolarization effect, all the domains in the horizontal grid edges align along the x-axis and those in the vertical grid edges align along the y-axis. However, every edge-crossing region is occupied by local $90^{\circ}$ domain structure (two $90^{\circ}$ walls).

It seems that for the grid FE lattice, the domain walls favor the self-adaptive evolution so that the electric polarizations neighboring the voids can align in parallel to the boundary rather than any other choice. This is a common feature almost applicable for all the cases. It is found that the domain splitting sequence can be realized by properly assigning the void array periodicity $(\lambda)$ and void size $L_{s}$ in a broad range, thus a flexible and easily accessed approach to control the domain structure is obtained. In comparison with FE islands, the grid FE lattice seems to offer rich possibilities of domain patterns [17].

In order to see more clearly the domain splitting, we choose the case of the $4\times 4$ void array and $L_{s}=6$ for an illustration. The snapshoted domain patterns at four annealing times $(T=0.05)$ are shown in figs. 2(a)–(d). Starting from the initial state in fig. 2(a), each stripe-like domain with vertical polarization begins to split at positions near the voids and generates two branches, as shown in fig. 2(b). The splitting develops continuously with time along the $\theta\sim135^{\circ}$ direction and ends near the neighboring void. This sequence thus generates the final well-aligned and regularly split domain structure shown in fig. 2(d) and fig. 1(d).

Zoom In Zoom Out Reset image size

Considering the well-organized split domain structures, we also consult the Fourier transform spectrum of the dipole correlation function in the wave vector $(k=(k_{x},\ k_{y}))$ space. The Fourier spectrum of the spatial correlation function is defined as

$$\begin{equation} S(k)=\frac{1}{L^{4}}\sum\limits_{ij} {P_{i} \cdot P_{j} \exp [ik\cdot (r_{i} -r_{j} )]}, \end{equation} \tag{ 10 }$$

where subscripts i and j denote the lattice sites. Figures 2(e)–(g) show, respectively, the Fourier spectra of the FE lattice free of void (fig. 1(a)), the FE lattices with $4\times 4$ void array and $L_{s}=6$ without any annealing and with sufficient annealing. A comparison of figs. 2(a) and (b) helps us to identify the indications of the spots. The bright spot at $k=(4, 4)$, as indicated by the arrow in fig. 2(e), reflects the stripe domain pattern whereas the weaker spot at $k\sim(12, 12)$ comes from the high-order effect from the Fourier transform. The additional dim spots in fig. 2(f) obviously come from the $4\times 4$ void array of $L_{s}=6$, noting that the brightest spot is at $k=(4, 4)$ from the stripe domain pattern. Figure 2(g) plots the Fourier spectrum of the well-developed domain pattern after the sufficient annealing. The brightest spot at $k=(8, 8)$ characterizes the well-split stripe domain pattern. This spot originates from the shifting of the brightest spot in fig. 2(f), due to the splitting of the stripe domains.

To understand the physics underlying such domain splitting behaviors, one may evaluate the spatial contours of the dipole-dipole interaction $(f_{\textit{dip}})$ and electromechanical energy $(f_{\textit{ela}}+f_{es})$. For convenience of comparison, figs. 3(a) and (b) show the spatial contours of $f_{\textit{dip}}$ around one void of $L_{s}=6$ at annealing time t = 0 and $t=5 \times 10^{6}\ \text{mcs}$, respectively. It is clearly shown that those polarizations perpendicular to the void surfaces at t = 0 re-align themselves to be parallel to the surfaces during the annealing, resulting in a significant reduction of the local $f_{\textit{dip}}$, which is one of the driving forces for the domain splitting. This process can also be illustrated by the spatial contours of the internal electric field. Figures 3(c) and (d) present the spatial contours of the magnitude of the internal electric field $\vert \boldsymbol{E}_{1}\vert$ at t = 0 and $5\times 10^{6}\ \text{mcs}$, with the arrows indicating the distributions of $\boldsymbol{E}_{1}$. At the up and left boundaries of voids, the dipole-dipole fields are initially offset by the strong depolarization field. The total internal field in this region is relatively small, as shown in fig. 3(c). The annealing process drives the dipoles to realign in parallel to the voids boundaries, leading to the modulated internal electric field, as shown in fig. 3(d) in comparison with fig. 3(c), focusing on the void-field boundaries. Figures 3(e) and (f) present the spatial contours of $(f_{\textit{ela}}+f_{es})$. Similar behaviors are observed and the existence of voids releases remarkably the local electromechanical energy. In particular, the high $(f_{\textit{ela}}+f_{es})$ energy on those domain walls (red color) in fig. 3(e) is fully released, and the whole region has much lower $(f_{\textit{ela}}+f_{es})$ energy via the polarization re-alignment. It is noted that the total gain in the electromechanical energy and dipole-dipole interaction is much greater than the cost for additional domain walls arising from the domain splitting. Therefore, the domain splitting is thermodynamically spontaneous. It should be mentioned that those energy terms in eq. (1) other than $f_{\textit{dip}}$ and $(f_{\textit{ela}}+f_{es})$ may change too upon the domain splitting. However, these changes are much less remarkable than $f_{\textit{dip}}$ and $(f_{\textit{ela}}+f_{es})$.

Zoom In Zoom Out Reset image size

One of the favored consequences of the domain splitting is the expected enhancement of the electromechanical performance. To confirm this prediction, we evaluate the piezoelectric coefficient $d_{33}$ defined as $d_{33}=\text{d} \eta^{[10]}/\text{d} E^{[10]}$ where $\eta^{[10]}= \langle\eta_{xx}\rangle=\langle Q_{11}P_{x}^{2}+ Q_{12}P_{y}^{2}\rangle$, $E^{[10]}$ denotes the electric field along the x-axis, and $\langle \rangle$ represents the configuration average [28]. When the averaging is made over the whole lattice including the voids, the evaluated piezoelectric coefficient is denoted by $d^{\textit{eff}}_{33}$, while the coefficient is denoted by $d_{33}$ if the averaging is made over the lattice excluding the voids. We take $-0.2<E^{[10]}<0.2$, within which term $\eta^{[10]}$ exhibits good linear dependence on $E^{[10]}$. The evaluated $d_{33}$ and $d^{\textit{eff}}_{33}$ as a function of $L_{s}$, respectively, are plotted in fig. 4(a). One observes that both $d_{33}$ and $d^{\textit{eff}}_{33}$ increase rapidly with increasing $L_{s}$ and reach the maximal as $L_{s}=4\text{-}8$. Afterwards, both $d_{33}$ and $d^{\textit{eff}}_{33}$ begin to fall and become low as $L_{s}$ is over 12. It is noted that more than 20 times enhancement in $d_{33}$ is obtained as $L_{s}$ increases from zero to $L_{s}=4$. This implies that the introduction of the void array allows a remarkably enhanced piezoelectric coefficient, a surprising but highly favored effect.

Zoom In Zoom Out Reset image size

To find the correlation of this piezoelectric enhancement with the domain structure, one plots the spatial contours of the local $d_{33}$ for the lattices of $L_{s}=0$ and $L_{s}=8$ in figs. 4(b) and (c), respectively, where the shadowed stripes roughly mark the domain walls. As shown in figs. 3(a) and (c), the stripe-like $90^{\circ}$ domain structure free of any void has both the high dipole-dipole interaction $(f_{\textit{dip}})$ and the high electromechanical energy $(f_{\textit{ela}}+f_{es})$. In particular, the high $(f_{\textit{ela}}+f_{es})$ implies that the whole lattice is in the seriously strained state. Therefore, additional lattice deformation becomes highly resistive, i.e. $d_{33}$ is very small. As expected, the $d_{33}$ contour for $L_{s}=0$ is roughly homogeneous over the whole lattice with small magnitude, and only the domain walls exhibit slightly higher $d_{33}$. Differently from this case, the lattice with the $L_{s}=8$ shows very high local $d_{33}$ not only on the domain walls, as illustrated in fig. 4(c). The physics is simple. On the one hand, the sufficiently released dipole-dipole interaction and electromechanical energy by the voids introduction benefit substantially to the high mobility and deformation capability of the domains under electric field. On the other hand, the domain splitting allows the doubling of the domain walls, contributing to the electromechanical response.

It should be mentioned that the present simulations are far from perfect. We only deal with the effects of the square void array on the $90^{\circ}$ domain structure and electromechanical behaviors. Also, only the 2D thin film lattice with in-plane FE domains is considered here. One may consider various arrays of voids of different shapes and deal with more realistic domain structures in FE thin films, where in-plane domains and out-of-plane domains coexist. The simulations may be extended to 3D cases with much more complicated domain structures. Nevertheless, for these extended cases, the underlying physics would be similar and the core issue here is how to efficiently release and even modulate the concerned free energy terms by implanting various geometric structures, so that the favored FE and electromechanical performances can be enhanced. In this sense, the present work, in addition to the FE islands, disks, dots, or composite structures, may shed light on alternative approaches to engineer the domain structure of FE materials for specific purposes. In fact, a recent experiment revealed that introduction of low permittivity defects into $\text{KTiOPO}_{4}$ lamellae can dramatically modulate the domain switching dynamics [29]. An anisotropic electric field pattern near the air holes is believed to be responsible for the observed results. In the present work, the internal electric field also shows a spatial distribution due to the depolarization field near the film-void interface. This experiment together with our findings here allows us to believe that modification of the FE domain structure by introducing voids-like defects is effective.

In summary, we have investigated the domain patterns and piezoelectricity of the tetragonal FE thin films with square grid structure. It is found that by properly choosing the size and array of the square voids, the stripe-like domains near the voids may spilt regularly. The competition between the dipole-dipole interaction and the electromechanical energy is responsible for the domain splitting. The calculations on the longitudinal piezoelectric coefficient predict a remarkable enhancement of the piezoelectric coefficient, caused by the sufficiently released dipole-dipole interaction and electromechanical energy due to the voids introduction and thus the doubling of domain walls.

This work was supported by the National Key Projects for Basic Researches of China (2015CB654602), the Natural Science Foundation of China (51431006, 51332006, 11374147), and the Priority Academic Program Development of Jiangsu Higher Education Institutions, China.