Electromechanical control of polarization vortex ordering in an interacting ferroelectric-dielectric composite dimer

TL;DR

This study uses a computational model to explore how a pair of ferroelectric nanospheres embedded in a dielectric medium respond to electric fields, revealing hysteretic and tunable polarization vortex behaviors influenced by particle separation and field orientation.

Contribution

It introduces a free-energy based model to analyze polarization vortex control in ferroelectric-dielectric composites, highlighting tunable responses via particle arrangement and external field direction.

Findings

System exhibits hysteretic and tunable responses to electric fields.

Polarization vortex configurations depend on particle separation and field orientation.

Potential for designing tunable composite metamaterials.

Abstract

Using a free-energy based computational model, we have investigated the response of a system comprising two interacting ferroelectric nanospheres, embedded in a dielectric medium, to a static external electric field. The system response is hysteretic and tunable by changing the inter-particle distance and the orientation of the applied field, which strongly modulates the field-driven long-range elastic interactions between the particles that propagate through the dielectric matrix. At small separations, the sensitivity of the system behavior with respect to the electric field direction originates from drastically different configurations of the local vortex-like polarization states in ferroelectric particles. This suggests new routes for the design of composite metamaterials whose dielectric properties can be controlled and tuned by selecting the mutual arrangement of their ferroelectric components.

Equations

$$\frac{\partial P_k}{\partial t}=-\Gamma \frac{\delta \mathcal{F}}{\delta P_k},$$

$$\frac{\partial }{\partial x_j}{C}_{ijkl}\left({\epsilon }_{kl}-Q_{klmn}{P}_m{P}_n\right)=0,$$

$$\frac{\partial }{\partial x_j}\left({ϵ}_b\frac{\partial \Phi }{\partial x_j}\right)=\frac{\partial P_k}{\partial x_k},$$

Full text

††preprint: This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. ††preprint: This article appeared in Appl. Phys. Lett. 113, 092901 (2018) and may be found at https://aip.scitation.org/doi/10.1063/1.5046080.

Ferroelectric-dielectric composite systems display a wide range of functional behavior. For particulate dispersed phases, precise control of the inter-particle spacing and arrangement Feng et al. (2017), particle shape morphology Hao et al. (2014); Caruntu et al. (2015), surface crystallography Wang et al. (2018) and concentration Hao et al. (2017), as well as selection of materials for both phases (with the matrix being of ceramicAirimioaei et al. (2017); Nenasheva et al. (2017) or polymeric originHuang et al. (2010); Paniagua et al. (2014); Hao et al. (2017); Feng et al. (2017)), can lead to designs with preprogrammed properties and novel functionalities. In particular, remarkable behavioral features, such as superparaelectricity Huang et al. (2010), enhanced dielectric tunability Zhou et al. (2008) and energy storage density Huang et al. (2010); Zhang et al. (2014); Liu et al. (2017), multi-state ferroelectric (FE) switching Gruverman et al. (2008); Li et al. (2017), and reduced conductivity and loss Airimioaei et al. (2017); Nenasheva et al. (2017) were already demonstrated in ferroelectric-dielectric composites, which makes these systems attractive for numerous future technological applications Alitalo and Tretyakova (2009).

A number of synthetic techniques have been recently used to fabricate nanoparticles of binaryJana et al. (2004); Jun et al. (2006); Tang and Gao (2016) and ternaryPark et al. (2005); Tyson et al. (2014); Papaefthymiou et al. (2015); Caruntu et al. (2015); Zhang et al. (2017) oxides in a variety of shapes and sizes. The properties of these particles, as well as the nature of FE ordering in them, have been characterized with multiple experimental approaches, including Bragg coherent diffraction imagingKarpov et al. (2017) and second harmonic generation microscopyCherifi-Hertel et al. (2017); Ma et al. (2017); Timpu et al. (2017) for sizes of above 50–60 nm, and aberration corrected transmission electron microscopy (TEM) Polking et al. (2012); Yadav et al. (2016) and local piezoforce microscopy for particles of to 5–15 nm in size Rodriguez et al. (2009); Polking et al. (2012). Both monodomain Polking et al. (2012) and vortex-like Rodriguez et al. (2009); Karpov et al. (2017); Li et al. (2017) polarization topologies were observed in the FE particles, with the latter patterns expected to be utilized in memory elements for new paradigms in computing Scott (2007); Gruverman et al. (2008).

Formation of topological polarization patterns in isolated FE particles has also been studied with different coarse-grained theoretical approaches, such as phase-field Wu et al. (2014); Liu et al. (2017); Mangeri et al. (2017); Yadav et al. (2016), semi-analytic Akdogan and Safari (2007); Glinchuk et al. (2008); Martelli et al. (2015); Qiu et al. (2015), and effective hamiltonian Naumov et al. (2004); Naumov and Bratkovsky (2008) methods. For example, in previous work Mangeri et al. (2017); Pitike et al. (2018), we evaluated size and shape dependent phase transitions in individual FE nanoparticles embedded in a dielectric matrix, observing a monodomain to vortex-like to polydomain sequence of transformations for the increasing particle size. However, in a composite system, this situation corresponds to a highly dispersed arrangement of particles within the matrix, with negligible electrostatic and elastic interactions between them. On the other hand, it is expected that more promising functional behavior can be obtained in complex systems consisting of many interacting particles, the aggregate responses of which should be heavily influenced by the particle size, shape, anisotropy and mutual arrangement Feng et al. (2017); Cai et al. (2017).

As accurate treatment of composite materials including numerous interacting FE particles presents a formidable computational challenge even on the coarse-grained level of theory, the next step towards developing a better understanding of such compounds involves investigation of systems with few interacting particles to contrast their behavior against the single-particle case. Here, using the time-dependent Landau-Ginzburg-Devonshire (TDLGD) approach Li et al. (2001); Mangeri et al. (2017), we explore a dimer system, consisting of two PbTiO3 (PTO) nanospheres embedded in a linear dielectric medium represented by SrTiO3 (STO). Our main goals include determination of a zero-field polarization pattern within the dimer as a function of the inter-particle separation $\eta$, followed by evaluation of the quasi-static response of the dimer system to an electric field that can be applied in different directions (as shown in the top panel of Fig. 1). Throughout this study, both particle diameters are kept fixed at $d=10$ nm: the size identified in our previous investigationMangeri et al. (2017) as the one that supports a zero electric field vortex-like polarization state.

Following the TDLGD approach, the Cartesian components $k=1,2,3$ of the local polarization field, $\mathbf{P}(\mathbf{r})$, are evolved in time as

[TABLE]

where ${\mathcal{F}}$ is the total electric enthalpy of the composite, and $\Gamma$ is a coefficient related to the FE domain wall mobility Hlinka (2007). Detailed expressions for all of the energy terms included in ${\mathcal{F}}$ are provided in the Supplemental Materials (SM), together with the utilized materials parameters for PTO and STO.

We assume that during the system time evolution, elastic strains relax at much faster time-scales than $\mathbf{P}$ Li et al. (2001, 2002). Therefore, at each time step, for fixed polarization, the strain fields $\varepsilon_{ij}=\frac{1}{2}\left(\partial u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i}\right)$ can be obtained from the mechanical equilibrium equation

[TABLE]

where $C_{ijkl}$ and $Q_{ijkl}$ are the elastic stiffness and electrostrictive tensor components, respectively, and $\mathbf{u}(\mathbf{r})$ is the elastic displacement field. The electrostatic potential $\Phi(\mathbf{r})$ is computed from the Poisson equation

[TABLE]

where $\epsilon_{b}$ is the background dielectric tensor Hlinka and Marton (2006), assumed to be isotropic.

We solve Eqs. (1–3) using finite-element methods within the open-source MOOSE frameworkGaston et al. (2009) and the Ferret module Fer developed by the authors.

In the absence of an applied external electric field, the boundary conditions on the six sides of the system volume (chosen to be sufficiently large with the dimer at its center) are $\mathbf{u}=0$, and $\Phi=0$. In the presence of an applied $\mathbf{E}$ field, a potential difference $\Delta\Phi=\Phi_{1}-\Phi_{2}$ is applied to two opposing sides of the box. The initial condition on $\mathbf{P}(\mathbf{r})$ is chosen as a random field of small local polarizations such that $\langle\mathbf{P}\rangle\approx 0$, i.e., mimicking the high-temperature paralectric state. As ${\mathbf{P}}$ is evolved in time [Eq. (1)], energy is dissipated from the system, with the simulation terminated when the relative change in ${\mathcal{F}}$ between consecutive time steps is approximately $0.01\%$.

An example of a dimer equilibrium state for $\eta=10$ nm under the $\mathbf{u},\Phi\to 0$ boundary conditions is presented in Fig. 1. The polarization field aligns tangentially to the particle surface to minimize surface charges, creating a vortex-like pattern in each of the two spheres. Such an equilibrium structure is different from the flux-closure points typically observed at intersections of four $90^{\circ}$ FE domain wallsMcQuaid et al. (2014); Wu et al. (2014). The pattern of Fig. 1 has diffuse domain walls, with the $\mathbf{P}(\mathbf{r})$ distribution continuously rotating about a central axis Mangeri et al. (2017); Karpov et al. (2017) represented by an axial director $\mathbf{v}$. The cylindrical vortex-core region penetrates through the sphere and usually has polarization magnitudes that are lower than at the surfaces. This behavior is different from that of ferromagnetic vortices, in which the magnetization density magnitude is constantLee et al. (2014); Zhou et al. (2015) well below $T_{\rm C}$.

For inter-particle spacings \eta\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}10 nm, the core-region orientations in the two particles can become strongly correlated. Specifically, the angles between the directors $\mathbf{v}$ of the two spheres can adopt a discrete set of values, suggesting that certain mutual core-region orientations within the dimer are more energetically preferred than others. In contrast with the single particle case, in some instances we also observe the development of large polarization, comparable to the surface values, in the vortex-core regions of both spheres. Some examples of such polarization patterns are shown in the SM for $\eta=2$ nm. For large inter-particle distances \eta\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}10 nm, we observe no pronounced correlations between the core-region orientations of the dimer particles; the orientation of the director v in each sphere depends mostly on the specific configuration of $\mathbf{P}(\mathbf{r})$ generated by the random initial condition. Therefore, in the absence of an externally applied $\mathbf{E}$, the dimer particles do not strongly interact with each other: the stray fringing electric fields from the surface and volume charge distributions of each sphere are effectively screened by the high dielectric permittivity ($\epsilon_{b}=300\epsilon_{0}$) of the surrounding matrix, and the strain fields propagating from the surface of each sphere into the matrix decay within a short distance.

We next turn to the analysis of the quasi-static dimer $P$-$E$ curves Mangeri et al. (2017); Boardman et al. (2005). In the previous study, the $P$-$E$ loops of an isolated PTO sphere with $d=10$ nm embedded in STO were found to exhibit shelves of roughly constant $\overline{P_{E}}=\overline{\textbf{P}\cdot\hat{\textbf{E}}}$ at smaller applied fields ($\overline{X}$ denotes volume-average of a scalar field $X$), which was associated with the presence of a vortex-like phase that transformed into a uniform monodomain phase when the field was increased, as shown in the insert of Fig. 2(b). A useful measure of the system vorticity is the FE-volume-averaged winding number density $\xi=|\textbf{P}\cdot\left(\nabla\times\textbf{P}\right)|$; $\xi\approx 0$ in the monodomain phase, and saturates to a constant value in the polydomain phase, however, in the adjoining vortex-like phase region it goes through a sharp maximum Mangeri et al. (2017).

The $P$-$E$ loops presented in Fig. 2 (a) and (b) are computed for $\eta=2,6,10$, and $22$ nm for both orientations of the applied electric field shown in the top panel of Fig. 1. For each point of the loop, the components of local $P_{E}(\mathbf{r})$ are averaged over both spheres. The averaged winding number density $\bar{\xi}$ is presented in panels (c) and (d) of Fig. 2 weighted with the average polarization along the field $\overline{P_{E}}$ as a function of the applied field magnitude $E$.

The orientation of the applied field relative to the dimer axis profoundly influences the shape of the $P$-$E$ response. With the field perpendicular to the axis connecting the particles ($\hat{\mathbf{E}}||\hat{\mathbf{z}}$) [Fig. reffig:hyst(a) and (c)], both the $\overline{P_{E}}$ and $\bar{\xi}\cdot\overline{P_{E}}$ dependencies exhibit only minor changes with increasing inter-particle spacing $\eta$. In particular, the transition between monodomain and vortex-like states occurs at the same value of $E\approx$ 0.2–0.3 MV/cm for all values of $\eta$. As the applied field switches sign, a small domain of correlated $\mathbf{P}$ along the direction of $\mathbf{E}$ begins to nucleate within the vortex-like phase in each sphere. With the increasing strength of the applied field, these domains grow until reaching saturation, which results in a gradual decrease and then complete disappearance of $\bar{\xi}$. This behavior seems to be universal and not strongly dependent on the value of inter-particle spacing $\eta$.

When the field is applied along the axis connecting the particles ($\hat{\mathbf{E}}||\hat{\mathbf{x}}$), as shown in panels (b) and (d) of Fig. 2, a radically different system behavior is observed, as manifested by the appearance of shelf-like features in both $\overline{P_{E}}(E)$ and $\bar{\xi}(E)$ curves, indicating the presence of the vortex-like phase. The transition into the vortex-like state happens at the same value of $E$ as in the previous case, i.e., 0.2–0.3 MV/cm. The width of the shelf-like feature on the $\overline{P_{E}}(E)$ curve is roughly inversely proportional to the value of $\eta$, as can be seen from the comparison of the loops for $\eta=2$ nm and $\eta=6$ nm [see Fig. 2 (b)]. The $\bar{\xi}(E)$ dependence exhibits sharp peaks on both ends of the shelf that are especially pronounced for $\eta=2$ nm and $\eta=6$ nm curves and can serve as good indicators for an onset of a phase transition. Inspection of the corresponding $\mathbf{P}(\mathbf{r})$ distributions shows that these peaks arise during a formation of multiple vortex-like patterns in each sphere. These patterns then merge into a single vortex with a core direction aligned with the direction of the field; some images illustrating the process are presented in the SM. For \eta\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}20 nm, the system $P$-$E$ loops become similar to that of an isolated sphere of identical size, cf. the $\eta=22$ nm curve and the curve in the inset of Fig. 2(b).

In order to better understand the origins of the varying dimer system response under different orientations of the applied field we examine the behavior of elastic strain fields as a function of the direction of $\mathbf{E}$. In Fig. 3, color maps of the diagonal components of the elastic strain tensor $\varepsilon_{ii}$ along the field are presented for the case of inter-particle separation of 6 nm and the field magnitude of 1.0 MV/cm applied in the $-\hat{\mathbf{z}}$ direction (left panel) and $-\hat{\mathbf{x}}$ direction (right panel). At this electric field strength, both particles in each configuration are in the monodomain state, with the polarization parallel to the applied field. Due to the nearly constant polarization, the elastic strain inside each sphere along the field direction $i$ is tensile, $\varepsilon_{ii}\approx 0.0165\%$. Therefore, the matrix region in between the particles experiences tensile strain for $\hat{\mathbf{E}}||\hat{\mathbf{z}}$ and compressive strain for $\hat{\mathbf{E}}||\hat{\mathbf{x}}$, as shown, respectively, in the left and right panels of Fig. 3. Elastic distortions extend for $\sim$10 nm from the sphere surfaces into the matrix and decay to zero at the boundaries of the computational domain (see SM for an example). The strain transfer mechanism in the composite arises from the electromechanical coupling between the FE particles and the dielectric medium: the applied electric field elongates the particles slightly along the direction of the field, which in turn causes elastic distortions within the matrix.

For the $\hat{\mathbf{E}}||\hat{\mathbf{x}}$ case, the strain in the region between the spheres is approximately additive for small separations $\eta$, which leads to a strong elastic coupling between the spheres. In turn, this results in a substantial decrease of $\overline{P_{E}}$, compared to the $\hat{\mathbf{E}}||\hat{\mathbf{z}}$ case: e.g., $\overline{P_{E}}\equiv\overline{P_{z}}=0.35$ vs. $\overline{P_{E}}\equiv\overline{P_{x}}=0.20$ C/$\mathrm{m}^{2}$ for $\eta=6$ nm and $E=0.7$ MV/cm, as shown in panels (a) and (b), respectively, of Fig. 2. The origin of this difference stems from contrasting orientations of the particle vortex cores with respect to $\mathbf{E}$, as shown in Fig. 4. In the $\hat{\mathbf{E}}||\hat{\mathbf{z}}$ case, shown in panel (a), core orientation vectors $\mathbf{v}$ of both spheres are perpendicular to $\mathbf{E}$ and it is the large local polarization from the surface regions that contributes to the average value of $P_{E}$. In contrast, in the $\hat{\mathbf{E}}||\hat{\mathbf{x}}$ case, shown in panel (b), the directors $\mathbf{v}$ of both spheres are aligned along E, with only the core-region polarization, which is greatly suppressed compared to the surface regions contributing to $\overline{P_{E}}$. Therefore, it is the $\mathbf{v}||\hat{\mathbf{E}}||\hat{\mathbf{x}}$ system orientation that gives rise to the shelf-like features with depressed near-constant $\overline{P_{E}}$ in Fig. 2(b) at small inter-particle separations. As the electric field increases, the core-region polarization saturates and then proliferates in the $y$-$z$ plane until a monodomain state is formed and vorticity vanishes.

We have investigated the influence dimer size and electric field orientation on the ordering and response of a FE-dielectric composite system. Our work demonstrates that the dimer size and electric field direction can be used to control the dielectric response in such composite systems, with potential applications in electronic and elecctromechanical devices.

Supplemental Material

See supplemental material for additional details on the calculations and supporting information.

Acknowledgments

JM acknowledges funding under the Ministry of Education, Youth and Sports of the Czech Republic, grant number CZ.02.2.69/0.0/0.0/16$\_$027/0008215. OH was funded by the Department of Energy, Office of Science Basic Energy Sciences Division of Materials Science and Engineering. JM is grateful to A. Jokisaari for assistance with calculations.