Phase Diagram of a Strained Ferroelectric Nanowire

: Ferroelectric materials manifest unique dielectric, ferroelastic, and piezoelectric properties. A targeted design of ferroelectrics at the nanoscale is not only of fundamental appeal but holds the highest potential for applications. Compared to two-dimensional nanostructures such as thin ﬁlms and superlattices, one-dimensional ferroelectric nanowires are investigated to a much lesser extent. Here, we reveal a variety of the topological polarization states, particularly the vortex and helical chiral phases, in loaded ferroelectric nanowires, which enable us to complete the strain–temperature phase diagram of the one-dimensional ferroelectrics. These phases are of prime importance for optoelectronics and quantum communication technologies. Author Contributions: Conceptualization, I.A.L.; methodology, I.A.L., V.M.V., S.K., M.S., Y.A.T. and L.B.; software, F.D.R., M.S., M.A.P. and Y.A.T.; validation, A.G.R., A.S., S.K. and V.M.V.; investigation, L.B., M.A.P., F.D.R., S.K. and A.S.; writing—original draft preparation, V.M.V., S.K., I.A.L. and A.G.R.; writing—review and editing, V.M.V., S.K. and I.A.L.; project administration, I.A.L., V.M.V., M.S. and


Introduction
Ferroelectrics are being turned into a leading material in modern electronics [1,2] owing to their exceedingly broad range of applications enabling innovative devices, such as memory and logical units, tunnel junctions, field-effect transistors, non-volatile ferroelectric random access memory, nano-electromechanical systems, energy-harvesting devices, advanced sensors and photocatalysis systems. A remarkable feature of ferroelectrics is that, in the course of the occurring dimensional downscaling of devices, they not only enable the superdensely packed circuits but also come up with new functionalities related to the emergence of topological polarization states and specific surface-induced electrodynamics [3][4][5]. The two-dimensional and quasi-two-dimensional ferroelectric thin films and superlattices harboring a variety of topological excitations such as soft domains [6], vortices [7], skyrmions [8], and merons [9] are now understood reasonably well. At the same time, cognizance surrounding the properties of topological excitations in their lowerdimensional relatives, zero-dimensional nanodots, is much less extensive, although similar excitations-vortices [10][11][12][13][14], skyrmions [15,16], and hopfions [17]-have been predicted. Furthermore, while the fabrication methods of one-dimensional nanorods, nanotubes, and nanowires have achieved considerable progress [18][19][20][21][22][23], the topological dynamics in these one-dimensional nanostructures remain even more puzzling.
The past progress in this study of ferroelectric one-dimensional nanostructures has been mostly dealing with general issues of the thermodynamic properties of the uniform polarization states [24][25][26][27]. The importance of "new physics" related to the emergence of unconventional topological vortex states was also realized [28], but only a few studies considered the systematics of their appearance [29][30][31]. In particular, the important practical aspects of operations with the different topological states remain underexplored.
Here, we investigated the emergent topological phase states in a strained ferroelectric nanowire arising due to the electrostatic confinement of polarization and demonstrate the ways of manipulating and switching these states by the load applied to the ends of the nanowire.
Here, we advance the existing Ginzburg-Landau-Devonshire (GLD) approach, via creating a synergistic technique harnessing the analytical treatment, the phase-field numerical calculations, and the atomistic simulations. A simultaneous use of these components not only enables new revelations but, most importantly, brings in strong reliability to derived results. These results show how new topological phases form in a confined strained nanoscale wire. We consider a nanowire made of PbTiO 3 (PTO), taking the latter as an exemplary ferroelectric material. We construct the strain-temperature phase diagram of the system and reveal the three qualitatively different polar phases: the vortex state (v phase) with polarization swirling around the nanowire c axis; the helical state (h phase), with the polarization screwing along the c axis; and the uniform polarization state extending along the c axis (c phase). Each state has unique physical properties for applications. In particular, the h phase hosts the spontaneously broken polarization chirality, easily tunable by the applied strain and temperature. This property of tunable chirality is of prime importance for optoelectronics and quantum communication technologies and does not have natural analogs.

Functional
A general approach to the description of ferroelectrics is based on the Ginzburg-Landau-Devonshire (GLD) functional, F = FdV, which in the case of a free-standing sample is the Gibbs thermodynamic potential of the system. The corresponding free-energy density, F(P i , σ ij , ∂ i ϕ), harnesses the ferroelectric, elastic, and electrostatic components, written in terms of the polarization vector, P i , stress tensor, σ ij , and electric field E i = −∂ i ϕ (where ϕ is an electrostatic potential), respectively. For the pseudo-cubic perovskite crystal, free energy density takes the form [32]: where the first brackets correspond to the uniform polarization energy [32], written in a form as in [33]; the next term accounts for the energy, related to the polarization gradient [34]; the stress-polarization coupling is described by the third term [32]; the fourth term represents the elastic energy contribution; and the coupling of polarization with the electric field and the electrostatic energy are given by the last two terms [35]. Indices i, j, k, l take the values 1, 2, 3 (or x, y, z), and the summation over repeated indices is assumed. For PTO with pseudo-cubic symmetry, the coefficients in free energy density (1) take standard values [36] (with cubic symmetry permutations) , a σ 11 = −0.73 × 10 8 C −4 m 6 N, a σ 12 = 7.5 × 10 8 C −4 m 6 N, a 111 = 2.6 × 10 8 C −6 m 10 N, a 112 = 6.1 × 10 8 C −6 m 10 N, and a 123 = −37 × 10 8 C −6 m 10 N. The components of the stress-polarization coupling tensor, Q ijkl , are Q 1111 = 0.089 C −2 m 4 , Q 1122 = −0.026 C −2 m 4 , and Q 1212 = 0.03375 C −2 m 4 . The components of the elastic compliance tensor, s 1111 = 8 × 10 −12 m −2 N −1 , s 1122 = −2.5 × 10 −12 m −2 N −1 , and s 1212 = 9 × 10 −12 m −2 N −1 are calculated as the inverse of the elastic tensor, s ijkl = c −1 ijkl . The gradient energy coefficients are G 1111 = 2.77 × 10 −10 C −2 m 4 N, G 1122 = 0, and G 1212 = 1.38 × 10 −10 C −2 m 4 N [34]. The constant ε 0 is the vacuum permittivity, and ε b 10 [37] is the background dielectric constant due to nonpolar ions.
The free energy density defined by Equation (1) written in terms of elastic stresses is ordinarily used for the description of the bulk free-standing systems. When discussing nanosystems confined in certain directions, it becomes more suitable to use the displacements corresponding to the confined degrees of freedom as driving parameters rather than stress variables. To that end, we go over to the strains, as to variables conjugated to the stresses, u ij = −∂F(P i , σ ij , ∂ i ϕ)/∂σ ij , and perform the Legendre transformation to replace part of the stresses by corresponding strains. As a result, we arrive at the mixed description in terms of both stress and strain variables, the choice depending on the particular geometry of the system and the property to be described.

Phase-Field Simulations
Phase-field simulations of the PTO nanowires are performed with the help of the FEniCS software package [38]. The evolution of the polarization vector field P i to its equilibrium fashion is governed by the relaxation equation: where the energy density is taken as a function of electric potential ϕ and the strain tensor u ij , obtained from the stress-dependent energy given by Equation (1) by the Legendre transformation [39]. Equation (2) is closed by the Poisson equation which takes into account electrostatic effects: and by the linear equation of elasticity: The value of the parameter γ is taken equal to one, as it is not crucial for the static case. The approximation of the time derivative on the left hand side of Equation (2) is accomplished by BDF2 variable-time stepper [40]. The system of nonlinear partial differential equations given by Equation (2) is solved by Newton's method. Linear Equations (3) and (4) are solved using the iterative generalized minimum residual method (GMRES) [41,42]. The computational space is represented by the cylindrical volume, partitioned into tetrahedrons by the 3D mesh generator gmsh [43]. All variables P i , ϕ, and u ij are periodically constrained in the z direction in addition to the assumed quasi-periodicity of the displacement along the cylinder axis, u z | z=0 = u z | z=l + hu 0 . Here, h is the cylinder height and u 0 is the applied strain. The procedure of calculation of the phase diagram is as follows. Let each point (T,u 0 ) of the system be initially at the paraelectric phase with a random distribution of P i with |P i | = 10 −6 C/m 2 components. Other initial states are the c phase with uniformly oriented polarization along the cylinder axis and the vortex v phase in which polarization rotates circularly around the cylinder axis. Then, let the system evolve with time and then choose for each point the state with the minimum energy defined by Equation (1).

Atomistic Simulations
Complementing the phenomenological approach, atomic-level descriptions are useful to model the microscopic behavior of ferroelectrics. Models based on interatomic potentials are computationally efficient for the simulation of finite-temperature properties in systems with a large number of atoms. In particular, a shell model approach with parameters fitted to first-principle results has been extensively used to study ferroelectric oxides [44][45][46]. In the shell model, atoms are thought to consist of an ion core coupled to an "electronic" shell in order to include its electronic polarization. The core and shell of an atom interact with each other as well as with the cores and shells of other ions via electrostatic interactions. The description also includes short-range interactions which are normally restricted to act between shells only [47]. The potential energy, V, is then a function of the positions of cores and shell, and is expressed as where V CS represents the contribution of the interatomic core-shell coupling, V LR contains the long-range electrostatic interactions among cores and shells of different atoms, and V SR accounts for the contributions of short-range interactions between shells. For the PTO model used here, cores and shells interact through an anharmonic spring V(w) = 1 2 k 2 w 2 + 1 2 k 4 w 4 , where w is the core-shell separation within an ion. Two different types of short-range potentials are considered: a Rydberg potential, V(r) = (a + br) exp(−r/ρ) is used for the Pb-Ti, Pb-O and Ti-O pairs, and a Buckingham potential, V(r) = a exp(−r/ρ) + c/r 6 , is used for O-O interactions. Coefficients k 2 , k 4 , a, b, c and ρ are the shell model fitting parameters. The multi-scale model with 23 independent parameters was originally developed to study bulk properties of PTO [48], and it was then successfully used to describe the behavior of surfaces, interfaces, and nanoparticles without including any additional coefficients [11,49,50].
The PTO wires are constructed from elementary perovskite unit-cells centered on Ti atoms with axes along pseudocubic directions (x, y and z axes lie along the [100], [010], and [001] directions, respectively). In addition, we consider that all surface faces consist of PbO planes. Molecular dynamics simulations with the atomistic model are carried out using the DL-POLY package [51] in a box with periodic boundary conditions. The nanowires extend along the z axis and they are surrounded by enough empty space to avoid any possible interaction between periodic replicas. Runs are performed by employing a constant volume and temperature (NVT) algorithm. The time step is set to 0.2 fs, which provides enough accuracy for the integration of the shell coordinates. The runs are made at temperature intervals of 50 • C. Each MD run consists of at least 100,000 time steps for data collection after 50,000 time steps for thermalization.
In contrast to the GLD functional, the atomistic description does not contain polarization as an explicit variable, and its value needs to be estimated from the positions of cores and shells. The local polarization is defined as the dipole moment per unit volume of a perovskite cell centered at the Ti atom, where v is the volume of the cell, z i and r i denote the charge and the position of the i-th particle, respectively, and w i is a weight factor equal to the number of cells to which the particle belongs.

Results
As an exemplary system, we employed the ferroelectric nanowire of the radius R = 4.8 nm and observed topologically different polarization textures for different strains and temperatures. Figure 1 summarizes the observed results as a u 0 -T phase diagram that exhibits the domains of existence of every discovered phase. The driving strain, u 0 = u zz , is calculated as the elongation of the lattice parameter of the nanowire c with respect to the relaxed bulk value c 0 as u 0 = c/c 0 − 1. The pairs of the cross-cut cylinders around the diagram exhibit the polarization structure in the observed phases. The left cylinder of each pair displays the polarization texture obtained from the phase field simulations visualized as the vector field, while the right cylinders show the results of the atomistic simulations visualized as polarization field lines. The red and blue lines in the phase diagram show the phase boundaries found from the phase-field and atomistic simulations, respectively. The high-temperature paraelectric phase is marked as the p phase in the phase diagram.
The vortex v phase arising at large compressing strains, harbors the closed field lines turning around the nanowire axis and lying in the plane perpendicular to this axis. The uniform c phase, arising at large tensile strains is the homogeneous phase when polarization is directed along the nanowire axis. Finally, the helical h phase is a new phase revealed by the phase-field simulations and arising when polarization exercises screw rotation around the nanowire axis. This phase exists at small strains and mediates the continuous transformation from a to c phases by changing the incline of the polarization vector with respect to the nanowire axis. The total view of the nanowire with the h phase is shown in Figure 1 on the left. In atomistic simulations, the h phase looks less regular, exhibiting a sort of propensity to history-dependent chaotic behavior. This may result from the surface or bulk-lattice pinning of polarization leading to the formation of the extra variety of the metastable states. Notably, as follows from the phase-field simulations, the second-order v-h and h-c transition lines join at high-temperatures, to form the line of the direct v-c firstorder transition along the short segment hitting the critical temperature line. Unfortunately, the precision of the atomistic simulation remains insufficient to test the thermodynamics in the region of the contact of the v, h, and c phases close to the critical temperature. Polarization, (C m ) P The low-right panel in Figure 1 demonstrates the practical way of the phase identification by plotting the average polarizations along and perpendicular the nanowire axis, P z = P 2 z 1/2 , andP ⊥ = P 2 x + P 2 y 1/2 as function of u 0 at given T. The regions wherē P z = 0,P ⊥ = 0,P z = 0,P ⊥ = 0 andP z = 0,P ⊥ = 0 correspond to the v, h, and c phase, respectively. This method gives satisfactory precision for the phase-field simulations (red curves) but is of more qualitative character in case the atomistic simulation (blue curves). As follows from Figure 1, the phase diagram obtained from the atomistic simulation is shifted with respect to the phase-field phase diagram to lower temperatures, and the corresponding characteristic values of polarization are smaller. The transition temperatures and polarization produced by the atomistic simulations are indeed underestimated with respect to those resulting from the phase-field simulations (and experimental results). This discrepancy is a consequence of the smaller equilibrium volume of the model and is related to the LDA database used to fit model parameters. It is well known that the volume and volume-dependent properties are underestimated when using the LDA approximation for the exchange-correlation functional in first-principle calculations which, in turn, lead to the significant underestimation of the polarization and transition temperatures. In bulk, for instance, the model gives a tetragonal ground state with the lattice parameter a = 3.859 Å, a tetragonal distortion c/a = 1.043, and the spontaneous polarization P = 0.54 C m −2 , while the experimental data are a = 3.90 Å, c/a = 1.065, and P = 0.75 C m −2 , respectively. The underestimation of the static structural properties is translated via the adjusted model to the finite temperature behavior. Molecular dynamics simulations showed the cubic-tetragonal transition at T c = 180 • C, which is of approximately 300 • C lower than the experimental value. Nevertheless, the qualitative temperature behavior of lattice parameters and polarization is correctly reproduced [48].

Discussion
Now we are in a position to incorporate our results into a general context of the existing research addressing the emergence of new polarization phases, homogeneous and inhomogeneous ones, in nanostructured ferroelectrics. A glorious avenue of exploring the properties of 2D and quasi-2D nanoscale ferroelectric heterostructures was paved by the seminal work of Pertsev, Zembilgotov, and Tagantsev [52], who constructed the phenomenological thermodynamic theory of strained ferroelectric films with uniform polarization and derived the strain-temperature phase diagram of the system. The applied mechanical in-plane strains were considered as the constraints for the partial Legendre transformations. Although the approach of [52] suffered from oversimplifications due to neglecting the depolarization effects and instabilities related to the emergence of ferroelectric and ferroelastic domains, it provided a gold standard for the initial approach to the consideration of the numerous nonuniform states discovered in ferroelectric 2D structures: complex networks of interplay of domains with 90 • and 180 • domain walls [53], periodic arrays of stripe soft polarization domains [35], having the texture of oppositely rotating vortices [7], or the double-periodic structure of cylindrical domains [54], having the textures of bubble skyrmions [8].
Here, we adapt the approach of [52] to 1D systems to gain, first of all, a general insight into the appearance of topological states in ferroelectric nanowires. We assume that the load applied to the ends of the nanowire results in the fixed strain u zz = u 0 along the nanowire (so that u 0 > 0 for stretching, and u 0 < 0 for compression). We further assume the absence of the off-diagonal strains, which implies u xz = u yz = u xy = 0. At the same time, the nanowire is supposed to be free to expand in the transversal x and y directions, while still being subject to the surface-applied pressure, p. This external pressure originates from either contact with the embedding material or from the Laplace pressure, p = µ/R, induced by the surface tension µ in the case of the cylinder sample [24,25]. Both cases result in the uniform in-plane compression stress σ xx = σ yy = −p.
The application of the partial Legendre transformations with the described constraints gives the effective free energy density F as Equation (6), that is similar to the thermodynamic potential obtained in [52] but with the differently renormalized coefficients denoted by asterisks: F = a * 1 P 2 1 + P 2 2 + a * 3 P 2 3 + a * 11 P 4 1 + P 4 2 + a * 33 P 4 3 + a * 13 P 2 1 P 2 3 + P 2 2 P 2 3 + a * 12 P 2 1 P 2 2 + a 111 P 6 1 + P 6 2 + P 6 3 + a 112 [P 4 1 P 2 2 + P 2 3 + P 4 2 P 2 1 + P 2 3 + P 4 3 (P 2 1 + P 2 2 )] + a 123 P 2 1 P 2 2 P 2 3 . (6) The analytical expressions for the renormalized coefficients and their numerical values for PTO are given in Table 1. Figure 2a presents the strain-temperature phase diagram (depicted by black lines) of the strained nanowires obtained as a result of minimization of the effective free energy density (6). The phase diagram contains an a phase with the transversal to the wire polarization, P = (P, 0, 0) (or equivalently the b phase with P = (0, P, 0)) and the c phase in parallel to the wire polarization, P = (0, 0, P); we assume here that the external pressure is absent, p = 0. Naturally, on the phase diagram, the negative (compressive) strain favors the a phase whereas the positive (tensile) strain favors the c phase.

Coefficient
Analytical Numerical, PTO Units Similar to the strained film [52], the behavior of the transition temperature of the ferroelectric nanowire, T c (u 0 , p) is depicted by the curve combining T ca and T cc dependencies on the applied strain, T c (u 0 , p) = max [T ca (u 0 , p), T cc (u 0 , p)], having a cusp at u 0 = 0 (at p = 0), see Figure 2a. Here, T ca and T cc are the Curie temperatures of the a and c phases, respectively; their dependence on u 0 and p, provided by conditions a * 1 = 0 and a * 3 = 0, is given in Table 2. The specific feature of the ferroelectric phase transition to the a phase is that it is the transition of the weak first order, because the fourth-order coefficients a * 11 , a * 22 are negative but very small in magnitude. We neglect the small discontinuity of this transition here. In the nonlinear region, below T c , the a phase occupies the substantial part of the diagram, including the clamped nanowire region with u 0 = 0. The c phase appears mostly at the large tensile strains and at high temperatures. The transition between the a and c phases occurs in a discontinuous way.
The preceding consideration, however, gives only a general structure of the phase diagram since it does not account for the electrostatic energy of the depolarization field induced by the surface bound charges arising at the surface points of polarization termination. Accordingly, the aforementioned phase diagram of the uniform ferroelectric phases is valid in the absence of depolarization effects, for instance, when the nanowire is submerged in the conducting media and the free conducting carriers screen the depolarization charges. Accounting for the electrostatic effects in the absence of screening modifies the emergent phases, making them nonuniform. The breakdown of the uniformity is caused by the necessity to reduce the value of the depolarization energy, ideally to zero, hence to avoid appearance of bound charges inducing the depolarization fields. Mathematically, this requirement implies that the system always prefers polarization textures with divergenceless polarization lines tangent to the lateral nanowire surface. The observed vortex v phase is one of such realized possibilities. The comparison of the uniform phase diagram with the results of our phase field simulations (shown by solid red lines in Figure 1, also replicated in Figure 2a by the red dashed lines) demonstrates that the nonuniform v phase forms from the uniform a phase, when the polarization vectors transversal to the wire swirl into the vortex to avoid the formation of bound charges at the lateral surface. The gradient energy, however, makes the v phase less favorable than the uniform a phase in case of screening. This reduces the critical temperature of the v phase and shifts its phase boundaries to the smaller strains with respect to the c phase which is evidently not affected by depolarization effects.
Importantly, the critical temperature of the vortex phase, T cv , can be analytically calculated by linearizing the GLD equations, obtained from the variation of the functional (1). The corresponding linearized equation in cylindrical coordinates r, θ for the axisymmetric polarization distribution P = P θ (r)θ assumes the form: which, taking the variational free boundary conditions ∂ r P θ (R) = 0, defines the critical temperature T cv and yields the corresponding vortex solution below T cv [29] as where J 1 is the first-order Bessel function, λ 1 = 1.8412 is the first zero of the derivative ∂ x J 1 (x), and ξ 0 = (G 1212 /AT ca ) 1/2 0.7 nm is the coherence length. The amplitude coefficient C is found by solving the nonlinear equations. Taking into account that on the u 0 -T diagram, according to Equation (8), the critical temperature of the v phase shifts down with respect to the screened a phase. We obtain that the cusp in the strain dependence of the total critical temperature of the system, T c = max(T cv , T cc ) also shifts down and to the left, to lower strains, with respect to the case of the uniform (screened) system. Shown by the solid red line in Figure 2a, this dependence perfectly coincides with the results of the phase-field simulations.
Another new phase, emerging from the strained ferroelectric nanowire is the helical h phase in which the polarization is inclined with respect to the nanowire axis and rotates around this axis. This phase, appearing at small strains, extends deep into the nonlinear region lying between the c and v phases and mediates the continuous transition between them by the gradual inclining of the rotating polarization vector from the direction parallel to the nanowire to the transversal to the nanowire direction. This helical phase is a state analogous to other topological states with non-zero helicity harboring skyrmions [15,16], coreless vortices [56], merons [9], and hopfions [17].
The importance of this phase is that it realizes a polarization state possessing a chiral asymmetry, a fundamental property of matter describing the systems distinguishable from their mirror images. The dependence of the quantitative characteristic of polarization chirality, the chirality density χ = P · (∇ × P) [17], on the strain and temperature is shown in Figure 2b,c, respectively, by the red lines. It is non-zero for the h phase and zero for the c and v phases, clearly indicating the region of the existence of the helical phase.

Conclusions
The theoretical investigation of a one-dimensional strained PbTiO 3 nano-wire undertaken herein allowed to reveal and model the emerging topological states in the ferroelectric as functions of strain and temperature. Two simulation techniques have been exercised and compared: the phase-field simulations employing the GLD functional and the atomistic simulations. Compressive and tensile strains were simulated at the upper and lower boundaries of the nano-wire along the axis. The strain-temperature phase diagram showing the different phases was produced and compared to the phase diagram that came from the simple free energy minimization procedure and which appeared more similar to the bulk phase diagram. Furthermore, the proposed novel methods enabled the discovery of two novel phases related to the considered depolarization effects, the vortex phase, having the polarization perpendicular to the wire axis, and the h phase, which mediates between the v phase and the uniform and previously known c phase.
As we demonstrated, the nanoscale confinement of a ferroelectric forming the strained nanowire stabilizes emerging topological states brings functionalities that do not exist in bulk materials. Nanowires supporting vortices and helical phase are of special interest for applications since these topological excitations are well controlled and manipulated by electric fields. This may promote and trigger giant piezo-and ferroelectric responses. Remarkably, we revealed that ferroelectric topological excitations in nanowires are controllable not only by electric fields, but are well tunable by the temperature and lattice strain. The possibility of switching between these states paves the way for using ferroelectric nanowires as base units for the multi-level logic ferroelectric devices [2,33,57], and neuromorphic computing circuits [1].
On top of that, chirality, a new functionality of confined ferroelectrics, has recently received explosive interest [16,17,58,59]. The newly discovered polarization chirality of nanowires, including its unprecedented tunability by the strains and temperature, opens a broad spectrum of fascinating directions promising to revolutionize different technological fields ranging from plasmonic and photonics relevant for quantum optoelectronic communications to chemical, pharmaceutical, and biomedical applications.