Phase-field simulations of vortex chirality manipulation in ferroelectric thin films

The ferroelectric chiral vortex domains are highly desirable for the application of data storage devices with low-energy consumption and high-density integration. However, the controllable switching of vortex chirality remains a challenge in the current ferroelectric community. Utilizing phase-field simulations, we investigate the vortex domain evolution and chirality formation in BiFeO3 thin films. By applying local surface charge or electric field, we demonstrate that the vorticity and the polarity can be manipulated by the initial bi-domain arrangement and the external field with different directions, respectively. By exchanging the domain arrangements, the opposite chirality can be obtained. Importantly, the topological vortex domain is retained after removing the external field. The vortex chirality can be switched reversibly with high reproducibility, which is beneficial to fatigue tolerance of the material in the operation. These results provide theoretical guidance for manipulating the vortex chirality in ferroelectric films.


INTRODUCTION
Chirality geometrically refers to the left-handed (LH) and righthanded (RH) phenomena 1 , which allows materials with chiral properties to be widely used, such as polarization optics 2,3 , sensors 4 , and stereochemistry 5 . Recently, nanoscale textures in ferroelectric materials have also been demonstrated to exhibit chirality 6 . Chiral topological defects, such as magnetic skyrmions 7,8 , polar vortices, and polar skyrmions [9][10][11] , have been reported, and these phenomena have attracted tremendous interest in material science and condensed matter physics 12,13 . Moreover, the relative insensitivity of these small-sized textures to external disturbances host the ability for their use in highly robust next-generation information storage devices 14,15 . This, in turn, has increased interest in the controlled manipulation of vortex chirality in ferroelectric materials.
Here, BiFeO 3 (BFO) serves as one of the most prototypical ferroelectric materials 16 . It is shown that its inherent rhombohedral (R) structure determines the eight polarization orientations along the pseudo-cubic <111> directions, and the adjacent orientations also determine the diversity of its domain walls (DWs), including 180°, 109°, and 71°DWs 17 . Previously, the chirality at DWs has been speculated due to the depolarization field and stress/strain inhomogeneity 18,19 . In the experiment, DW chirality has been proved in triclinic ferroelectric crystals 20 . Then, chiral 71°a nd 109°DWs appear in BFO films owing to the different polarization rotation ways at DW, which is confirmed by using the reciprocal and real-space characterization techniques 21 .
Theoretical works have reported that curled electric field (E) controls vortex chirality by switching toroid moment to form a new vortex in ferroelectric nanoparticles 22,23 , and the interaction of polarization and strain promotes vortex formation. However, with the local radial E stimulation, ferroelectric vortex pairs were created or annihilated in films by selecting the appropriate polarization location [24][25][26] , and characteristic shaped or graded composition [27][28][29] . In addition, the mechanical field was proven to control vortex morphology [29][30][31] . As for the interface strain, it also promotes polarization rotation to form a pair of opposite vortices 32 . Among them, the influence of flexoelectricity on polarization switching and DW formation has been extensively explored 18,33,34 . These results indicate that vortex evolution is sensitive to flexoelectricity [35][36][37] , and the flexoelectric coupling coefficients can be quantified by phase-field methods, where the vortex morphology is reflected 38 .
Based on the phase-field simulation, this work focuses on the formation mechanism and controllable manipulation of chiral vortices in a single domain and bi-domain in BFO thin films under local surface charge or electric field. The simulation results demonstrate that the bi-domain arrangements determine the inplane vorticity. The chirality is related to the initial bi-domain arrangements. Importantly, the chirality is reversibly switched by applying alternating surface charge or electric field. The vortex chirality remains unchanged after removing the external field. This study advances our understanding of vortex domain evolution and chirality manipulation in ferroelectric thin films.

RESULTS AND DISCUSSION
The emergent chirality in the vortex The schematics of four types of chirality defined according to the out-of-plane polarity and in-plane vorticity of the polarization vectors (P) [39][40][41] are shown in Fig. 1. The out-of-plane polarity points downward, that is, P z < 0 (marked as a black cross), which is defined as OP = −1. When the out-of-plane polarity is upward, that is, P z > 0 (marked as double black circles), then is defined as OP = +1. In addition, counterclockwise vorticity is CCW = −1 and clockwise vorticity is CW= +1, and the black arrows indicate the projection of the polarization in the x-y plane. The chirality is then defined as the product of the polarity and the vorticity. Accordingly, LH chirality is defined as LH = 1, which includes a CCW (−1) vorticity with downward (−1) polarity, and a CW (+1) 1 vorticity with upward (+1) polarity, as shown in Fig. 1a, d, respectively. Similarly, RH chirality is defined as RH = −1, which includes a CW (+1) vorticity with downward polarity (−1), and a CCW (−1) vorticity with upward polarity (+1), as shown in Fig. 1b, c, respectively. Figure 1e is a schematic of the eight polarization variants in R BFO. According to the symmetry, there are four upward (R1+, R2+, R3+, and R4+) and four downward (R1−, R2−, R3−, and R4−) polarization variants, respectively.
The vortex chirality in the single domain First, we apply positive/negative surface charge on the top surface of the 10 nm BFO thin film with an initial single domain, and the vortex can be generated as shown in Fig. 2a, b. The results of the other four initial single domains (R1−, R2−, R3−, and R4−) are shown in Supplementary Fig. 1. Figure 2c shows that nine positive surface charge regions are applied on BFO thin film with an initial R1+ domain. Interestingly, the vortex chirality is random with the initial single domain due to the degenerate energy state of polarization. Then, it should be noticed that the downward polarization switching appears in the local region with applied surface charge, whereas the domains with upward P z appear around the vortex. Here, the applied local positive surface charge is equivalent to a downward built-in electric field. The positive surface charge favors the negative polarization bound charge, and the in-plane domains with outward center-divergent polarization directions are formed 13,25 . In addition, the out-of-plane polarization points downward due to the negative built-in electric field. Although we cannot controllably manipulate the vortex chirality in a single domain, the discrepancy in polarization switching of different initial domains can promote us to recognize the formation mechanism of vortex domains. The percentages of four negative R domains induced in the single domain are shown in Fig. 2f. The results show that the percentage of R domains formed by 71°polarization switching (e.g., R1+ → R3−) is larger than that formed by 109°(e.g., R1+ →R2−, R4−) and 180°(e.g., R1+ →R1−) switching. Here, the height of the energy barrier determines the difficulty of the path of polarization switching 42 .
Then, we further explore the influence of surface charge density (σ surf ) on the domain evolution, as shown in Fig. 3. The domain structures under applied different positive σ surf are shown in Fig. 3a, and the corresponding energy densities evolution is shown in Fig. 3c. R1+ single variants are initial states and the R domains with downward out-of-plane polarization (P z < 0) are induced by applying positive charges. When the σ surf is 0.08 Cm −2 , R3− domain is dominant. As the σ surf continues to increase, the R1− domain region gradually increases. When the σ surf is 0.8 Cm −2 , a vortex with evenly distributed 4 R domains forms, as shown in Fig. 3a. The square-shaped local positive surface charge region is applied on the top surface, and the downward domain switching occurs in the charged region. Then, the surrounding domains (P z > 0) are formed to minimize the stray electric field. The polarization switching mechanism affects the R domain ratio in a single vortex, and 71°ferroelastic switching is the easiest. The surface charge densities determine the domain switching path, which in turn affects the vortex domain morphology. At the same time, the vortex core moves toward the center of the film with the σ Ã surf increases, as shown in Fig. 3b. The corresponding energy density (f) evolution is shown in Fig. 3c, where the decrease in f electric is due to the contribution of the built-in electric field caused by the increase in surface charge.
The expression of f Landau (Eq. 3) is closely related to polarization, where an increase in surface charge induces more polarization switching and leads to the increase of DW density, thereby increasing the Landau energy. When 71°DW converges to form a vortex, the vortex core exhibits continuous polarization rotation. To minimize the free energy and electrostatic energy of the system, the domain configurations transform to a curved DW in a metastable state, as shown in Supplementary Fig. 2.

The vortex chirality manipulation in bi-domain
The DW is the key to breaking the radial symmetry of the external field, which is conducive to the formation of more interesting topological vortex patterns 24,43 . A bi-domain arrangements with initial R3+|R3− polarization states are separated by 180°DW oriented along the y axis, and LH chiral vortex with CCW rotation can be obtained in Fig. 4a, as indicated by the polarization vector plot. Then, by reversing the initial bidomain arrangement from R3+|R3− to R3−|R3+, an RH vortex with CW rotation is induced as shown in Fig. 4b. We also consider bi-domain arrangements with initial R3+|R2− states separated by a 109°DW oriented along the y axis, as shown in Fig. 4c, and the reversed domain arrangement with initial R2−|R3+ states is shown in Fig. 4d. The LH chiral vortex domain with CCW rotation can be obtained, and the polarization projection in the x-y plane is shown in vector plots, respectively. The polarization orientation of the initial bi-domain determines vortex chirality. In other Here, the inherent DW chirality makes contributions to the charge-induced vortex chirality. In R3+|R3−, the polarization components (P y ) in two adjacent domains have opposite directions and show CCW vorticity seen from the right to the left so that the charge-induced vortex is CCW. In R3−|R3+, P y shows CW vorticity so that the charge induces a CW vortex. The two P y are in the same direction in Fig. 4b, d, so the two vortices have similar rotation directions. However, the charged 109°DW (R3+|R2−) are mechanically incompatible, but the final states of these DWs are stable in the modeling due to enough screening charges in short-circuit electrostatic boundary conditions. To understand the effect of mechanical compatibility of initial DW on vortex chirality, the charge-induced chiral vortices are obtained in neutral and mechanical permissible 109°DWs (R1+ |R2− and R3 −|R4+) 44,45 , respectively, as shown in Supplementary Fig. 6. In R1 +|R2− and R3−|R4+, the opposite vortex chirality is obtained after exchanging the bi-domain arrangements owing to the opposite P y in adjacent domains. The inherent P y rotation direction determines the vortex chirality, these results further confirm the above discussions.
To study the influence of DW orientation on the vortex chirality, we switch the DW orientation along the x axis, as shown in Fig. 4e-h. Opposite chirality is obtained by reversing the bidomain arrangement from R3+/R3− to R3−/R3+, respectively, as shown in Fig. 4e, f. Similarly, reversing the initial R3+/R2− to R2−/R3+, RH and LH vortices are obtained, respectively, as shown in Fig. 4g, h. Similarly, a pair of controllable chiral vortex domains, including polarity and vorticity type, can also be successfully achieved in the initial tri-domain combination, as shown in Supplementary Fig. 7.
The chirality formation is closely related to the initial bi-domain arrangements, the detailed evolution is shown in Supplementary  Fig. 8. As we know, four R domains positions are fixed due to the crystal orientation in R-BFO films 46,47 . In R3+|R3-bi-domains [ Supplementary Fig. 8a], at left-panel, R1− domain is larger than R4−, by 71°and 109°switching. On the right panel, positive σ surf determines the larger R3− domain. Then, the straight four R domains gradually evolve into an LH vortex. However, an RH vortex is shown in R3−|R3+ bi-domains [ Supplementary Fig. 8b   vorticity. The charge-induced vortex chirality is closely related to the initial bi-domain.
Although the successful application of surface charge-induced polarization switching in the electrochemistry 48,49 , the common path to achieve domain switching is the electric field or the electric bias in the experiments 25,46,50 . Therefore, the local radial electric fields created by applying square-shaped surface potential are studied on the bi-domain BFO films with 180°DWs, the surface potential distributions are shown in Supplementary Fig. 9. The results are consistent with the case of surface charge, a specific vortex chirality can be obtained by applying a local external electric field in the ordered initial bi-domain.

The strain distribution in the chiral vortex
In ferroelectric materials, there is an electrostrictive coupling effect between ferroelectric polarization and intrinsic strain (ε 0 ij ). Therefore, the elastic strain (e ij ) distribution can reflect the ferroelectric polarization distribution to a certain extent (Seen in Eq. (7)). We calculate elastic strain (e ij ) in chiral vortex induced by initial R3+ |R3−, R3−|R3+, R3+|R2− and R2−|R3+ bi-domain, respectively. Figure 5e-h shows elastic strain distributions along the corresponding red dashed line of Fig. 5a-d, respectively. For all the strain distributions, the axial strains e 11 , e 22 , and e 33 change significantly near the vortex core, e 11 and e 22 are tensile strains and e 33 is compressive. By comparing Fig. 5e, f, the shear strains e 23 and e 13 show different evolution behavior. Therefore, we plot the shear strain (e 23 and e 13 ) distribution along the 2D x-y plane in Fig. 5i-l. The e 23 and e 13 have maximum and minimum values with the opposite signs on both sides of the core, and we add the arrows to indicate the direction of e 23 and e 13 . In Fig. 5i, j, both e 23 and e 13 all have opposite directions, so the vortices show CCW rotation and CW rotation, respectively. Similarly, in Fig. 5k, l, e 23 and e 13 have the same direction on both sides of the core, so the vortices are all CCW rotation. The shear strains (e 23 and e 13 ) have closely related to the vortex vorticity in the core.
The external field manipulation of vortex chirality To explore the influence of the flexoelectric effect, the bi-domain combinations of R3+|R3− of 180°DW are selected as the initial states, as shown in Supplementary Fig. 10. The vortex chirality remains LH by considering the flexoelectric effect. Then, the strain gradient of the in-plane elastic strain (e 11 , e 22 ) and out-of-plane elastic strain (e 33 ) along the z axis slightly increase due to the flexoelectric coupling coefficients. Here, the flexoelectric effect contributes to the morphology of the vortex domain but does not change the chirality of the ferroelectric vortex.
To further achieve chirality switching, we demonstrate the reversible control of vortex chirality by alternating external fields. The schematics of alternate application of σ surf and E are shown in Fig. 6a, b, respectively. The reversible chirality of the vortex domain can be obtained by applying ± σ surf and ±E, as shown in Fig. 6c, d, respectively. These results demonstrate that the out-ofplane polarity changes with the sign of the external field in the vortex, while the vorticity is unchanged. Here, the width of the charge region (W c ) is 80 nm in Fig. 6, The stray electric field induces the outer domain with opposite P z direction. The adjacent inner and outer domains are 71°out-of-plane polarization orientation, which decreases the large out-of-plane depolarization field to protect the inner vortex domain stability. Then, to study the vortex stability, the topological chiral vortices are relaxed after we remove the surface charge. The morphology of the vortex domain shrinks slightly, while the chirality of the vortex remains stable, the surrounding domains remain stable and the 71°out-ofplane orientation relationship is still in the inner vortex domain and surrounding outer domains, as shown in Supplementary  Fig. 11. To understand the effect of the surface charge region size on the stability of the chiral vortex, we apply two widths of the charge region (W c = 40 nm, 60 nm) on the 180°DWs, as shown in Supplementary Fig. 12. When we remove the surface charge, the chiral vortex in a small charge region is unstable and the vortex merges with the surrounding domains and becomes a striped domain. Here, the outer large domains make the inner small vortex merge and disappear. Therefore, the vortex stability is closely related to the large charge region size in BFO films.
In short, the ferroelectric vortices can be induced in initial single and bi-domain regions under local surface charge and electric field stimulation. In a single domain, the initial domains determine the percentage and morphology of the vortex domain, the vortex formation requires sufficient surface charge density. As the charge density increases, the vortex core moves to the center. Our results show that vortex chirality is controllably manipulated by selecting the initial bi-domain arrangement with typical DWs. The inherent DW chirality determines the charge-induced vortex chirality. In other words, if the initial bi-domain polarization component's direction along the DW are opposite, the opposite vortex chirality can be obtained. The sign of the external field determines the vortex polarity, thus successfully regulating chirality. After we remove the external field, the vortices are stable only in the region with previously applied large enough charge density. These results provide a basis for understanding the chirality manipulation and vortex evolution in ferroelectric thin films.

METHODS
In the phase-field simulation, the ferroelectric polarization evolutions are described by the time-dependent Ginzburg-Landau equation [51][52][53] : where P i is the component of the polarization along with the [1 0 0] (x axis; i = 1), [0 1 0] (y axis; i = 2), and [0 0 1] (z axis; i = 3) crystal orientations in Cartesian coordinates, r and t are the spatial position vector and the simulation time, respectively, L is related to DW mobility, and f total is the total free energy density in the thin film, which is determined by the Landau free energy density (f Landau ), gradient energy density (f gradient ), electrostatic energy density (f electric ), and elastic energy density (f elastic ), and flexoelectric energy density (f flexo ): We define f Landau as the following sixth-order polynomial, a 1 , a 11 , a 12 , a 111 , a 112 , and a 123 are the dielectric stiffness and high order stiffness coefficients 54 . According to the Curie-Wiess law, a 1 depends on the temperature, and it can be defined as a 1 = (T−T 0 )/(2ε 0 C 0 ), where T is temperature, T 0 is the Curie temperature, the vacuum dielectric constant is ε 0 = 8.85 × 10 −12 Fm −1 , and C 0 is the Curie constant.
The flexoelectric energy density f flexo can be expressed as where f ijkl represents flexoelectric coupling coefficients and e ij is elastic strain. For cubic symmetry, f ijkl can be simplified as f 11 , f 22 , f 44 in Voigt's notation. Here, the flexoelectric coupling coefficients are set as: f 11 = 2.5 V, f 22 = 2.5 V, f 44 = 0.05 V, based on reasonable estimation 55 . Then, the driving force from the flexoelectric effect (E flexo i ) can be calculated as follows: The elastic energy density f elastic can be written as c ijkl e ij e kl ¼ 1 2 c ijkl ðε ij À ε 0 ij Þðε kl À ε 0 kl Þ; ði; j; k; l ¼ 1; 2; 3Þ where c ijkl , e ij , ε ij , and ε 0 ij is the elastic stiffness tensor, elastic strain, total strain, and eigenstrain of BFO films, respectively. The eigenstrain is ε 0 ij ¼ Q ijkl P k P l , Q ijkl is an electrostrictive coefficient tensor and can be marked as Q 11 , Q 22 , and Q 44 due to the symmetry. Here, there are three elastic constants of c 11 , c 12 , and c 44 in Voigt's notation for a cubic symmetry. The top surface of films is stress-free and the bottom is constrained by the substrate. In the bottom surface of the film, the equal biaxial misfit strain is considered along the x and y directions. The elastic stress is σ ij ¼ c ijkl e ij ¼ c ijkl ðε ij À ε 0 ij Þ, and the mechanical equilibrium equation is Here, on the top surface of the film σ i3 j x3¼nf ¼ 0.
The electrostatic energy density f electric can be written as where ε 0 is the vacuum dielectric constant, and k ij represents the background dielectric constants 56 , and here the isotropic k ij is taken as k 11 ¼ k 22 ¼ k 33 ¼ 45, which is consistent with the adopted value of the previous references 57 . Electrostatic energy comes from two contributions, one comes from spontaneous polarization, and the other comes from energy storage in vacuum space. The contributions of polarization come from the ionic displacements and background dielectric constants. E i ¼ À∇ i φ ¼ À ∂φ ∂xi is an electrostatic field, φ is electric potential. The shortcircuit boundary conditions are applied to solve the electrostatic equilibrium equation, which avoids the depolarization effect and bound charge screening effect. Then, the surface charge density (σ surf ) is applied on the top surface of thin films and the free charge density has the following boundary condition: ρ jz¼nf ¼ σ surf , and ρ ¼ 0 in the thin film. Assuming that the σ surf is that the charge occupies each grid point, σ surf ¼ σ Ã surf electrons=S grid area , one grid size is S grid_area = 1.0 × 10 −18 m 2 and one electron is equal to 1.6 × 10 −19 C. Accordingly, the electrostatic potential φ(r) at the thin film surface can be solved by electrostatic equilibrium (Poisson) equation 58 , where n f is the thickness of the thin film. The details of boundary conditions can be seen in ref. 33 .
In the simulations, the three-dimensional sizes of the simulated systems are 64Δx × 64Δx × 16Δx for a single domain, where Δx = 1 nm, and 128Δx × 128Δx × 16Δx for bi-domains, and 240Δx × 240Δx × 16Δx for tri-domain. In addition, the thickness of BFO film is 10 nm (n f = 10Δx), and the material coefficients of BFO are listed in Table 1.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.