Phase-field simulations of topological conservation in multi-vortex induced by surface charge in BiFeO₃ thin films

1. Introduction

A variety of exotic topological defects have been discovered in ferroelectric materials including flux-closure quadrants,¹ vortexes,^2,3 skyrmions,⁴ Merons,⁵ center-type domains⁶ and so on, which contain a wealth of interesting physical phenomena and new functions.^7,8,9 Through precisely controlling these topological structures, high-efficiency and multifunctional nanoelectronic devices^10,11,12,13 can be designed. As one of the typical topological structures, polar vortex has been created in oxide superlattices, nanostructures and polymers^{13,14,15,16,17,18,19,20} and has been observed to vibrate at extremely high frequencies.^21,22 Besides, the vortex core exhibits an unusually stable state with low energy,^17,23 making vortex an ideal unit for ultrafast, dense and low-crosstalk data processing and information storage applications.^21,22,24,25

Polar vortex can be a spontaneous formation in plenty of ferroelectric materials.^{24,26,27,28,29} For example, Huang et al. discovered unique hierarchical vortexes within polycrystalline ferroelectric materials.³⁰ Besides, by using electric field,^{22,31,32,33,34,35,36} surface charge,^37,38 temperature stimulation,^11,31,39 mechanical stress,^{40,41,42,43,44} strain,^10,45,46,47 oxygen vacancy regulation^43,48 and flexoelectricity,^49,50,51,52 people have been able to create and manipulate polar vortex precisely. However, considering the applications in actual devices, the topological conservation and interaction between the adjacent vortexes need to be further investigated.

In this work, we first explore the topological conservation of the surface charge-induced vortex evolution by constructing a phase-field model^36,38,53 in BiFeO₃ (BFO)^54,55,56,57 thin films. Multiple vortexes are induced by local surface charges in BFO films with different sizes, and the conservation of topological vortex is analyzed by calculating the winding numbers. Subsequently, the minimum distance between the adjacent surface charge regions without domain interaction is explored, and finally the corresponding storage density of the vortex memory is calculated. Our work provides a theoretical basis for the practical application of high-density vortex memories.

2. Theoretical Methods

In phase-field simulations, the spatial distribution of the ferroelectric polarization vector $P_i$ is a natural order parameter, which is used to describe the phase transition and the domain structure evolution. The time-dependent Ginzburg–Landau (TDGL) formula (1) and the semi-implicit Fourier transform method^58,59 can be used to solve the polarization field :

where $P_i$ ($i=1,2,3$) is the polarization component severally along the x-axis [110], y-axis [010] and z-axis [001] in the Cartesian coordinate system, $r$ and t are the spatial position vector and evolution time of the system respectively, L is the dynamical coefficient related to the domain wall motion, $F$ is the total Free energy of the system including Landau energy $f_{\mathrm{\text{Landau}}}$, gradient energy $f_{\mathrm{\text{gradient}}}$, elastic energy $f_{\mathrm{\text{elastic}}}$ and static electric field energy $f_{\mathrm{\text{electric}}}^{}$.⁶⁰ The Landau energy $f_{\mathrm{\text{Landau}}}$ can be written in the form of a polynomial, where, ${\alpha }_{ij}$, ${\alpha }_{ijkl}$ and ${\alpha }_{ijklmn}$ are the coefficients of the Landau polynomial related to the dielectric stiffness of different orders, which can usually be expanded to the order four through eight. We expanded to sixth order in this paper. The gradient energy $f_{\mathrm{\text{gradient}}}$ can be described as follows : where $G_{ijkl}$ is the gradient coefficient, $P_{i,j}=\frac{\partial P_i}{\partial x_j}$. The elastic energy $f_{\mathrm{\text{elastic}}}$ can be described as follows: where $C_{ijkl}$ is the elastic coefficient tensor, $e_{ij}$ is the elastic strain, ${\epsilon }_{ij}$ (${\epsilon }_{ij}^0=Q_{ijkl}{P}_k{P}_l$) and ${\epsilon }_{ij}^0$ are the total elastic strain and intrinsic strain of the system, respectively. $Q_{ijkl}$ is the electrostrictive coefficient, which represents the coupling degree between polarization and lattice freedom in the ferroelectric. The electrostatic field energy $f_{\mathrm{\text{electric}}}$ in this system can be expressed as follows : Here, $E_i$ is the electrostatic field, $k_{ij}$ and ${\epsilon }_0$ are dielectric constant and permittivity of vacuum of the material, respectively. $E_i=-{\nabla }_i\phi =-\frac{\partial \phi }{\partial x_i}\phi$, $\phi$ is the potential. Short-circuit boundary conditions were adopted for the solution of the electrostatic equation, which could effectively avoid the depolarization effect and the shielding effect of the bound charges. In this work, in order to induce the polar vortex, a specific region of surface charge density (${\sigma }_{\mathrm{\text{surf}}})$ was applied to the top of the film. The application of surface charges to stabilize the topological domain structure had been studied and verified in previous experiments.⁶¹ The boundary condition is ${\rho }_{|z=nf}={\sigma }_{\mathrm{\text{surf}}}$, and $\rho =0$ in a thin film. ${\sigma }_{\mathrm{\text{surf}}}$ represents that charges occupy each grid point, ${\sigma }_{\mathrm{\text{surf}}}={\sigma }_{\mathrm{\text{surf}}}^{\ast }\times \frac{electrons}{S_{\mathrm{\text{gri}}{\mathrm{\text{d}}}_{\mathrm{\text{area}}}}}$, here each electron is charged according to $1.6\times {10}^{-19}\mathrm{\text{C}}$. And ${\sigma }_{\mathrm{\text{surf}}}^{\ast }=5$ is chosen because a uniformly distributed 4R vortex domain can be induced in this case.⁶² The size of the charge region is $50\mathrm{\Delta }x\times 50\mathrm{\Delta }y(\mathrm{\Delta }x=\mathrm{\Delta }y=1\mathrm{\text{nm}})$, the square of the grid is $S_{\mathrm{\text{gri}}{\mathrm{\text{d}}}_{\mathrm{\text{area}}}}=1.0\times {10}^{-18}{m}^2$.

In this work, the associated coefficients of BFO were in agreement with the simulation parameters used in Ref. 38 by Liu et al. (Table 1). A square surface charge region with sides of 50nm is placed on the top of the film to induce a large domain vortex as shown in Fig. 1(a). The three-dimensional simulation size is $128\mathrm{\Delta }x\times 128\mathrm{\Delta }y\times 20\mathrm{\Delta }z$ ($\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z\mathrm{\text{=1nm}})$. As four large vortexes are induced, four identical surface charge regions are placed uniformly in a larger size measuring $160\mathrm{\Delta }x\times 160\mathrm{\Delta }y\times 20\mathrm{\Delta }z(\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z=1\mathrm{\text{nm}})$ shown in Fig. 1(b). In addition, the thickness of the film was 10 nanometers ($n_f=10\mathrm{\Delta }x$) and the temperature was at 300K in our study.

3. Results and Discussion

To describe the topological characteristics of polar vortex, the winding numbers of the polarization vector in a local area are usually calculated.^12,63,64,65 When the circle is traversed counterclockwise in the parameter space of polarization, the angular increment in the counterclockwise direction is written as positive (otherwise negative). The winding number n is defined as the change of angle between the x axis and polarization vector along the whole loop divided by 2$\pi$, which requires the absolute value of the angle to be at least 2$\pi$ if n is nonzero. n equals 1 indicates that there is a vortex in the domain, or an anti-vortex is displayed if n${}_A$ equals −1.^12,66

To examine the conservation law of surface charge-induced vortex in BFO films, we investigated the vortex evolution in the films with initial random eight irregularly distributed R domain variants (see in Supplementary Information Fig. S1(a)) and initial R1+ domain variant (polarization along the [111] crystalline direction, see in Fig. S1(b)), respectively. For a more detailed analysis of the pattern, we plotted the domain morphology and the corresponding polarization profiles on the top and bottom layers of the film, respectively. We also calculated and plotted the distribution maps of winding numbers in different conditions through MATLAB program. In Fig. 2, vortex and anti-vortex are induced in BFO films with initial random domains using a positive surface charge region. In Fig. 2(a), pinwheel-shaped 4 R vortex domains surrounded by a cluttered array of other R domain variants (hereafter referred to as variant(s)) can be observed in the middle region on the top surface of the film. The downward polarization state is energetically favorable since the positive surface charge induces a built-in electric field along the $-z$ direction in Fig. 2(b). The yellow dot indicates vortex and the winding number is marked as n = 1, while the green dot indicates anti-vortex and the winding number is marked as n${}_A$=−1 in Fig. 2(c). There is an interesting phenomenon that the vortex at the center of the bottom surface of the film occupies a tiny area, which appears to be newly grown from the original vortex domain in Figs. 2(d)–2(f), possibly related to the short action time and long action distance from the surface charge. We can conclude that the number of vortex and anti-vortex on both the top and bottom surfaces is 1 as the initial domain is random. We also calculated the winding numbers of topological domains with initial random domain and R1+ domain layer by layer (see in Fig. S2), respectively. The value of winding numbers in each layer is 0. No vortex appears on any other surfaces of the topological domains except the top and bottom because of the size effect of the films.

As shown in Fig. 3, we continued to induce vortex domains using surface charge in BFO film with initial R1+ domain. Fig. 3(a) displays an elliptical vortex domain on the top surface. Here, the R2+ and R4+ variants grow around the R4- and R2- variants in the vortex domain, respectively. In Fig. 3(b), the polarization vector changes in a clockwise direction, consistent with the vortex rotation. Furthermore, the enlarged polarization view shows an anti-vortex in the lower left corner, which is the same position as the green dot in Fig. 3(c), so that the winding number is marked as n = 1, n${}_A=-1$. A smaller vortex is observed in the center of the domain on the bottom film like Fig. 2(d). There are a pair of vortex and anti-vortex observed clearly in Figs. 3(e) and 3(f), and the winding numbers can be marked as n = 1, n${}_A$ = −1. Hence, the vortex and anti-vortex will appear in pairs on the top and bottom surfaces of the film induced by a single surface charge region, and the sum of the winding numbers certainly should always be 0.

We further verify the topological conservation law when multi-vortex domains are induced simultaneously based on the above research results. Here, four distinct and neatly aligned vortex domains are induced in Figs. S3 and S4. The pinwheel domain consists of R1-, R2-, R3- and R4- domain variants, which were encircled by R3+, R4+, R1+ and R2+ variants, respectively to maintain the equilibrium. It’s worth mentioning that the formation of stable vortexes induced by surface charges in thin films with initial random domains requires a longer evolution time. Every time a vortex is successfully induced, there is also an anti-vortex in its lower right corner in Fig. S3 (the initial domain is random) or lower left corner meanwhile in Fig. S4 (the initial domain is R1+ single domain). This is the further evidence for the simultaneous existence of the vortex and anti-vortex during the introduction of positive surface charges.

We also found a very interesting pattern during the domain evolution by adjusting the $\mathrm{\Delta }l$ in turn from 20nm to 2nm between the surface charge regions. Four distinct disjointed vortexes were induced at the beginning. Then, R2+ and R4+ variants, respectively, surrounding the adjacent vortexes began to stick together. Afterwards, different R variants that made up the vortex in the central region of the film began to connect with each other. Finally, four big vortexes coalesced into one larger vortex while a tiny anti-vortex existed. The critical distance $\mathrm{\Delta }l$ that determines the vortexes interaction characteristics is different with different initial domains. To facilitate understanding, the morphologies of the vortexes in different cases are shown. The positions of (anti-) vortexes on the top surface of the thin films are located by calculating the winding numbers. The special evolution of the topological domains as the initial domain is random is shown in Fig. 4(a). The R2+ and R4+ variants surrounding the adjacent vortexes haven’t touched each other when $\mathrm{\Delta }l\ge 12$. The vortexes still exist independently at $\mathrm{\Delta }l=9$nm. However, the domain variants which formed the vortex begin to touch each other when $\mathrm{\Delta }l$ is 8nm. When $\mathrm{\Delta }l$ is shortened to 2nm, four vortexes completely overlap to form a larger vortex that covers four charge regions. Besides, three pairs of tiny vortexes and anti-vortexes coexist under the large vortex. The winding number is marked as $n=4$ and $n_A=-4$ in Fig. 4(b). It is worth mentioning that this is because of numerical iteration, tiny vortexes are generated and identified in the drawing, which doesn’t affect the vortex conservation.

However, the R2+ and R4+ variants enclosing the vortexes begin to adhesion to each other when $\mathrm{\Delta }l$ is shortened to 15nm in the thin film with the initial R1+ single domain in Fig. 5(a). The R3$-$ variants which formed the vortex overlap with each other in the center of the film when $\mathrm{\Delta }l=10$ nm in Fig. S4. And at $\mathrm{\Delta }l=2$nm, a larger square vortex is formed. The winding number in Fig. 5(b) began to change with the $\mathrm{\Delta }l$ is 8nm, which is marked as $n=5$ and $n_A=-5$. In the end, a pair of vortex and anti-vortex exist which marked as $n=1$ and $n_A=-\mathbf{1}$ and consistent with the topological morphology. Therefore, all vortexes will merge into one larger vortex and remain stable when the surface regions are close enough.

For more in-depth analysis of the evolution of vortexes with $\mathrm{\Delta }l$, the statistics of the winding numbers on the top surface of the thin films as the initial domains are random (Fig. 6(a)) and R1+ domains (Fig. 6(b)) respectively were calculated and plotted. The two numerical curves representing the winding numbers of vortexes and anti-vortexes range from 1 to 5 and are perfectly symmetric distributed, which means that the sum of the winding numbers is 0. As the distance $\mathrm{\Delta }l$ changes, the winding numbers also changes regularly. When the initial domain is random, it is intuitive to see the number of vortexes is more stable than the initial R1+ variant. The winding numbers remain $n=4$ and $n_A=-4$ respectively as the distance $\mathrm{\Delta }l$ is greater than 6nm when the initial domain is random, indicating high stability. However, the distance $\mathrm{\Delta }l$ to maintain a stable state is more than 8nm when the initial domain is R1+ variant.

We aim to investigate the conversation law of vortex domains induced by surface charges in BFO thin films and discuss the practical significance of designing BFO vortex memories. It can be seen that the induced vortex domain distribution is more stable when the initial domain is random, and the corresponding minimum distance $\mathrm{\Delta }l$ is 9nm from the above qualitative analysis of vortex domain morphology and winding numbers. In the real data storage device, it is well known to use the ratio between the storage of the node data itself and the total storage of the entire node architecture, i.e., the storage density, to highlight the storage performance. In BFO films, it is appropriate to use the surface density, which is the measurement of the number of information bits in the surface area, to characterize the storage density. The basic unit is bit⋅inch${}^{-2}$, where 1 bit can only store 0 or 1. Therefore, the surface density of the film, when expressed in GB⋅inch${}^{-2}$, can be calculated using formula (7). The storage density is at least 172.609GB⋅inch${}^{-2}$ when the vortex domain is induced by the square surface charge regions with sides of 50nm. As early as 2016, the size of the independent storage unit used to store data has been reduced to 10nm,⁶⁷ with the density requirement for memory reaching TB⋅inch${}^{-2}$. Although our simulation results in this paper haven’t met the contemporary demand, the small size of the domain itself allows the induction of even smaller domains through the reduction of the surface charge area to achieve higher density memories, which is our next step.

4. Conclusions

In summary, the topological conservation in multi-vortex induced by square surface charge regions of 50nm side length in BiFeO₃ thin films based on the phase-field simulations has been studied. By mapping the corresponding vortex domain morphologies, we found that the vortex domains are distributed like pinwheel-shaped islands. Meanwhile, in the evolutions of vortex domains induced by surface charges, it is proved that the vortex and anti-vortex exist in pairs and the total winding numbers is always 0. Through adjusting the distance $\mathrm{\Delta }l$ between adjacent surface charge regions, we also found that the induced (anti-) vortex distribution is relatively stable when $\mathrm{\Delta }l$ exceeds 8nm if the initial domain is random or $\mathrm{\Delta }l$ exceeds 10nm if the initial domain is R1+ variant. As the distance $\mathrm{\Delta }l$ is reduced to 2nm, four vortexes on the film converge into a larger tetragonal vortex with a tiny anti-vortex at the corner. We can further explore more complex and denser multiple vortex arrangements by adjusting the size of the films and the surface charge regions. This will provide a basic understanding and a reference point for the further exploration of high-density vortex memories.

Acknowledgments

This work was supported financially by the National Key Research and Development Program of China (Grant No. 2019YFA0307900), the National Natural Science Foundation of China (Grant Nos. 51972028 and 12004036) and the China Postdoctoral Science Foundation (Grant No. 2020M680375). This project is also supported by Young Elite Scientists Sponsorship Program by BAST (No. BYFSS2023072) and State Key Laboratory of New Ceramic and Fine Processing Tsinghua University (No. KFZD202201). The phase-simulation results in this work were obtained using the software package Mu-PRO (www.mupro.co).

ORCID

Mengyuan Wang  https://orcid.org/0009-0007-4660-9830

Di Liu  https://orcid.org/0009-0000-2718-0775

Jing Wang  https://orcid.org/0000-0002-9099-7339

Deshan Liang  https://orcid.org/0000-0002-1179-3902

Houbing Huang  https://orcid.org/0000-0002-8006-3495

Supplemental Materials

The Supplemental Materials are available at: https://www.worldscientific.com/doi/suppl/10.1142/S2010135X23430026.