Electric field driven evolution of topological domain structure in hexagonal manganites
K. L. Yang,1Y. Zhang,1S. H. Zheng,1L. Lin,1Z. B. Yan,1J.-M. Liu,1,2,3and S.-W. Cheong4
1Laboratory of Solid State Microstructures, Nanjing University, Nanjing 210093, China
2Institute for Advanced Materials, South China Normal University, Guangzhou 510006, China
3Hubei Normal University, Huangshi 435000, China
4Rutgers Center for Emergent Materials and Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854, USA
(Received 28 May 2017; published 9 October 2017)
Controlling and manipulating the topological state represents an important topic in condensed matters for both fundamental researches and applications. In this work, we focus on the evolution of a real-space topological domain structure in hexagonal manganites driven by electric field, using the analytical and numerical calculations based on the Ginzburg-Landau theory. It is revealed that the electric field drives a transition of the topological domain structure from the type-I pattern to the type-II one. In particular, it is identified that a high electric field can enforce the two antiphase-plus-ferroelectric (AP+FE) domain walls with=π/3 to approach each other and to merge into one domain wall with=2π/3 eventually if the electric field is sufficiently high, whereis the difference in the trimerization phase between two neighboring domains. Our simulations also reveal that the vortex cores of the topological structure can be disabled at a sufficiently high critical electric field by suppressing the structural trimerization therein, beyond which the vortex core region is replaced by a single ferroelectric domain without structural trimerization (Q=0). Our results provide a stimulating reference for understanding the manipulation of real-space topological domain structure in hexagonal manganites.
DOI:10.1103/PhysRevB.96.144103
I. INTRODUCTION
Topological states/defects and relevant emergent phenom- ena have been for a long period attracting attention in theoret- ical condensed-matter physics and mathematical physics, and recently become hot topics due to the finding of a set of novel topological quantum states described in momentum space [1–5]. A topological state/phase can be defined either in the momentum space or real space. Distinctly different from the topological physics in momentum space such as topological insulators and Weyl semimetals, real-space topological states of interests in condensed matters and materials have been relatively localized, and the local topological defects and related phase transitions (e.g., Kosterlitz-Thouless transitions) were discussed in the pioneer works of Kosterlitz and Thouless [6]. Nevertheless, the real-space topological physics has been much less discussed, and so far most of the earlier investigations focused on characterizations of topological defects/states including morphology, topological invariances, order parameters, and stability, in systems like liquid crystals, frustrated magnetic systems, graphene, and even carbon nan- otubes [7–12]. Symmetry consideration of these topological states has been also an often discussed issue [1].
In fact, possible functionality of real-space topological states/defects is yet a realm to be explored. So far, synthesis of an isolated topological defect can be a challenge and its detection and manipulation have been rarely investigated either. In the overall sense, the response of a localized topological state and its dynamics against either intrinsic or extrinsic stimuli/perturbations remains to be a topic of less interest. A well-known character of a topological state is its topology robustness measured by a topological invariance [1].
This invariance does not necessarily imply any robustness of geometric pattern which can be largely flexible against those stimuli, noting the misleading issue that a topological pattern
is indestructible due to its topology protection. Along this line, a major issue or opportunity is to search for possible topological states, no matter whether they are localized or extended, which can be given emergent functionalities of application potentials. Here we discuss a multiferroic material with real-space topological state and investigate the evolution of this state to external electric field.
It has been found that hexagonal manganites RMnO3 (R=Sc,Y,Dy-Lu), upon proper pretreatments, may exhibit a specific kind of real-space topological state, noting that hexag- onal manganites are well-attracted multiferroics [13–15].
A distinct difference but also an interesting point here is that the topological states stem from the intercoupled ferroelectric (FE) domains and antiphase (AP) domains, allowing the close connection between the materials’ functionality and the topological state. While each topological defect may be localized, a high density of such defects could self-organize into a large scale domain pattern or structure, here called a topological domain structure [16], consisting of vortex- antivortex pairs, as schematically shown in Fig. 1(a) for a guide of eyes. In this sense, an investigation on the response of such a topological structure against stimuli/perturbations becomes concerned.
For the details of topological structure in hexagonal man- ganites, one may consult relevant literature [17–24]. Briefly, each topological defect (vortex or antivortex) results from a trimerization-type structural phase transition induced by size mismatch between rare-earth ionR and Mn-O layers. Below certain temperature (Tc), this trimerization is triggered by periodic tilting of the MnO5trigonal bipyramids from the high- symmetryP63/mmc structure. Consequently, the generated lattice distortion makes the unit cell three times larger in size and the symmetry is lowered down to the polarP63cmspace group. For a straightforward understanding of such topological
FIG. 1. A schematic drawing of real-space topological structure in hexagonal manganites RMnO3, projected on theab plane: (a) type-I pattern and (b) type-II pattern. The structural trimerization leads to the threefold antiphase domain structure (α,β,γ), i.e., the Z3 symmetry. The ferroelectric polarizationP has two degenerate momentsP↑andP↓oriented respectively along the±caxis, with the Z2symmetry. The vortex-antivortex pairs are properly labeled.
structure, the phenomenological Landau theory treating this trimerization was proposed, and this theory starts from two order parameters: trimerization amplitude Q describing the tilting magnitude from the caxis and trimerization phase describing the tilting orientation [25]. This trimerization is disabled ifQ=0. Given the hexagonal lattice symmetry, three degenerate types of structural antiphase domains (α,β, andγ), described by theZ3symmetry, are available. Each antiphase can supply two different ferroelectric polarizations P (P↑ and P↓, indicating the upward and downward polarizations along the caxis), respectively corresponding to the upward distortion of two thirds of the Y ions and the downward distortion of one-third of Y ions or vice versa, set by the in or out orientation of the MnO5 trigonal bipyramid tilting.
In this way, the antiphase domains and ferroelectric domains are mutually interlocked, constituting the AP+FE domain structure. The two-state polarization (Z2symmetry) plus the three-state antiphase (Z3symmetry) allows the presence of six antiphase and ferroelectric domains (α+,α−,β+,β−,γ+,γ−) emerging from one point, constituting a topological defect called a vortex/antivortex described by the Z2×Z3 sym- metry. The vortex and antivortex correspond respectively to two degenerate winding orders, i.e., (α+,β−,γ+,α−,β+,γ−) and (α+,γ−,β+,α−,γ+,β−) [17,26]. In the projected two- dimensional plane (e.g.,abplane), each vortex must neighbor with one antivortex, constituting one vortex-antivortex pair, as shown in Fig.1(a)too. A number of such pairs form the complex topological pattern in real space, which has been discussed in the framework of graph theory [16]. Subsequently, various trials in order to track the pattern variants and their evolutions have been performed, suggesting that the domain
structure can be highly dependent of sample processing details including thermal treatments, imposed strains, and even surface oxygen deficiencies, etc. [13,27–30]. It was also reported that an electron beam is sufficient to nucleate an opposite domain from a parent domain [31]. These phenomena are more or less local and they may resist against or benefit the formation of topological states.
However, for a topological state, the topological robustness is an intrinsic property. For instance, it was reported that a shear stress could deform the vortex-antivortex pattern seriously into a stripelike pattern, owing to the so-called Magnus-type force which pulls the vortex-antivortex pairs apart in opposite directions and thus unfolds the network [28]. Nevertheless, these two types of domain structures have identical topological invariants. External stimuli do deform the pattern seriously but for most cases the topological invariance is preserved, although the functionality responses could be different upon the different patterns. What should be addressed here is that the topological pattern evolution does not show a simple scenario in response to various stimuli. One often observed sequence is the transition from the type-I domain pattern to type-II pattern.
Here the type-I pattern refers to a class of domain structures with roughly equal fractions ofP↑- andP↓ domains and the type-II pattern to those with one dominant polarization, as schematically shown in Fig. 1 for a guide to the eyes. The type-I pattern has roughly no net polarization but the type-II pattern exhibits a nonzero net polarization along one of the±c axes [16,29]. This transition was at least confirmed by electric field poling experiments but no details have been available [32].
We are interested in the topological structure in response to electric field along thecaxis without losing the generality.
This could be informative for realizing an electrocontrol of the topological structure. Due to the interlocking of FE domains and AP domains, we assigned the interlocked structure as the AP+FE domain. Two neighboring AP+FE domain walls with phase difference=π/3 could eventually fuse into one antiphase (AP) only wall with=2π/3, driven by electric field. Certainly, this fusing behavior should be highly resistive since the disappearance of aP↓ domain may probably be accompanied with a breakdown of the topological structure which has to overcome a extremely high energy barrier. Consequently, it would be useful if the evolution of the topological structure and its topological behavior can be tracked consecutively. In particular, the question of whether the topological structure can be suppressed by certain stimuli, such as electric field, is interesting.
In this work, we mainly address such an evolution including the dynamic process, using extensive phase-field simulations based on the Landau theory. First, we do “observe” the electric field driven fusion tendency of two neighboring AP+FE walls with =π/3 into one AP-only wall with =2π/3.
Second, the lattice details with the AP-only walls would benefit to our understanding of the topological structure evolution.
Third and surprisingly, sufficiently high electric field may trigger a transition of the lattice structure at the walls and vortex/antivortex cores from the low-symmetry trimerized lattice into nontrimerized lattice while they both may be ferroelectric. In particular, the AP+FE domain structure can be decoupled by sufficiently high electric field.
The remaining part of this paper is organized as follows.
In Sec. II, we discuss the variation of order parametersP, Q, and with varying electric field in the framework of Landau phenomenological theory. The method and procedure of phase-field simulations are briefly described. The main
results of the dynamic evolution of topological structure will be presented in Sec. III, and in particular how the vortex/antivortex core is destabilized by sufficiently high electric field is discussed. A brief conclusion will be given in Sec.IV.
II. PHENOMENOLOGICAL THEORY AND SIMULATION A. Landau theory
The Landau phenomenological theory in terms of the variations of trimerization and polarization of hexagonalRMnO3was developed recently and extended to explain a set of observed phenomena [25,33,34]. The topological state can also be properly predicted. The system free energy (density) is expressed
F =FL+FG, (1)
in which the Landau energyFLand gradient energyFGcan be written respectively as FL= a0
2 Q2+b0
4 Q4+c0
6Q6+c0
6 Q6cos 6−gQ3Pcos 3+g
2Q2P2+aP 2 P2, FG= 1
2
l=x,y,z
sQl
∂lQ∂lQ+Q2∂l∂l
+sPl∂lP ∂lP
, ∂l= ∂
∂l, (2)
where order parametersQandare the amplitude and phase of trimerization mode (K3), andP is the local polarization (or the local amplitude of polar mode 2−). The subscripts (l=x,y,z) are the major coordinate axes and (a0,b0,c0) are the constants for a free-energy polynomial onQextended up to the sixth order,c0is the anisotropic coupling factor betweenQand,gis the nonlinear coupling factor between the modeK3 and mode 2−,gis the coupling factor betweenQandP. It is noted that aP is the self-energy factor ofP extended up to the second order, while the case forP extended up to the fourth order will be discussed in Sec.III. Here coefficientsslQandslP scale the energy costs for the spatial variations ofQandP and they can be defined as the stiffness parameters forQandP respectively.
The simplest way to investigate the dependences ofQ,, andP onEis to add the electrostatic energy term to Eq. (1):
FE= −E·P = −(ExPx+EyPy+EzPz)⇐
E=Ex−→
i +Ey−→
j +Ez−→ k P =Px−→
i +Py−→ j +Pz−→
k
, (3)
where we only deal with polarizationP =Pzalong thec(z) axis, settingPx =Py =0. Certainly, it is sufficient to considerEz
in this case. HereP is scaled by the polar mode amplitude, giving the unit of ˚A, according to the treatment of Artyukhinet al.
[25], and the unit ofEzis eV/A, which is not the common unit of electric field. In order to obtain the real electric field for a˚ comparison with experimental data,Ez=EE0is set, whereE0=1.0 eV/A and˚ Eis a dimensionless number. The real electric fieldEr=2EE0/(9.031e)=22146Ewith the unit of kV/cm because an earlier analysis suggested that this unit of 1.0 ˚A for the polar mode amplitude is equivalent to an electric fieldE0=22 146 kV/cm for hexagonalRMnO3. Unfortunately, the value ofE0 is much larger than the measured coercive field,∼40 kV/cm for YMnO3as an example. This gap seems to be attributed to several reasons associated with the imperfections of samples under measurements although the existence of this gap does not change the physical essence of problems discussed here.
B. Ground states and effect of electric field
We first discuss the ground states under an electric field, without inclusion of domain-wall energy associated withFG. In this case, the system total energy is F =(FL+FE), and its minimization with respect to Q, , and P would result in the relationships among them. In proceeding, we obtain the ground states by solving the following set of equations:
∂F
∂Q =[a0+b0Q2+(c0+c0cos 6)Q4−3gQPcos 3+gP2]Q=0,
∂F
∂ =(−2c0Q6cos 3+3gQ3P) sin 3=0,
∂F
∂P = −gQ3cos 3+(gQ2+aP)P −Ez=0. (4)
A consequential algebraic formulation from Eq. (4) allows the order parameters to satisfy a0+b0Q2+(c0+c0)Q4−3gQcos 3gQ3cos 3+Ez
gQ2+aP
+g
gQ3cos 3+Ez gQ2+aP
2
=0,
=0,±π/3,±2π/3,π, P = gQ3cos 3+Ez
gQ2+aP
. (5)
Clearly, six degenerate sets of order parameters can be obtained from Eq. (5), which correspond to the six degenerate domain states. In more detail, one seesP ∝cos3atE=0, and thus the six states can be divided into two types: one type with upward polarization (P↑) and the other type with downward polarization (P↓). Correspondently, Eq. (5) can be rewritten into two sets:
a0+b0Q2+(c0+c0)Q4−3gQgQ3+Ez
gQ2+aP +g
gQ3+Ez
gQ2+aP
2
=0, =0,±2π/3, P = gQ3+Ez
gQ2+aP(P↑,P >0), a0+b0Q2+(c0+c0)Q4−3gQgQ3−Ez
gQ2+aP +g
gQ3−Ez gQ2+aP
2
=0, = ±π/3,π, P = −gQ3−Ez
gQ2+aP(P↓,P <0).
(6)
Now Eq. (6) can be used to calculateP as a function ofE, as shown in Fig.2(b)where twoP(E) branches are plotted, one forP >0(P↑) and the other forP <0(P↓). TheP↓state becomes unstable asEexceeds 1.11, and theP↑state becomes unstable asE < −1.11. This explains why the two curves pass across the P =0 axis and then terminate nearby. The calculated twoQ(E) branches are plotted in Fig.2(a)where the solid line parts represent the stable solutions and the dashed parts are the unstable solutions. Looking at the stable solutions, one finds that the variation ofQ(E) asP andEhave the same sign is slightly weaker than that asP andEhave the opposite signs. The evaluated energy density F(E) curves as P >0 andP <0 are plotted in Fig.2(c), where the data atQ=0 are inserted too for a comparison. In fact,Q=0 implies the disappearance of the six antiphase states, i.e., absence of the trimerization. The unstable solutions correspond to the regions in the sufficiently high field where theQ=0 state is lower than theQ =0 states in the total energy, to be discussed later.
Here several main features associated with the ground state can be highlighted. First, in the low-field range, amplitude Q is insensitive to E but the trimerization could lose its stability in the high-field range. Once such a destabilization occurs, the trimerization and thus the six antiphase states will be completely suppressed and replaced by the Q=0 state, but the ferroelectric state remains in this case. This suggests that the topological structure would be replaced by a normal
ferroelectric state in sufficiently high electric field. Second, a roughly linear dependence ofP onE, no matter forP >0 and P <0 states, is shown, suggesting no polarization saturation in the high-field range. This is certainly a prediction by the present Landau theory, probably inconsistent with realistic observations in some cases. Third, the ground state in the high-field range is no longer the antiphase structure ofQ =0, but replaced by theQ=0 state whose lattice structure remains elusive so far. This replacement is obviously a first-order phase transition, as evidenced by theQ(E) hysteresis in Fig.2(a).
C. Dynamic evolution of order parameters
Given the obtained ground state, one is now able to track the topological structure evolution using the phase-field simulations in which the polar coordinates Q and are transformed into the Cartesian coordinatesQx andQy where Qx =Qcos and Qy =Qsin [35], and the domain-wall energy FG is taken into account. Following the standard procedure in the literature [36,37], we start from the temporal evolution of order-parameter fields described by the time- dependent Ginzburg-Landau equations:
∂η(r,t)
∂t = −L δF
δη(r,t), η=Qx,Qy,P , (7) whereLis the kinetic coefficient. For details, one has
∂Qx(r,t)
∂t =L
sQx ∂2Qx(r,t)
∂x2 +sQy∂2Qx(r,t)
∂y2 +sQz ∂2Qx(r,t)
∂z2 − ∂fL
∂Qx(r,t) ,
∂Qy(r,t)
∂t =L
sQx ∂2Qy(r,t)
∂x2 +sQy∂2Qy(r,t)
∂y2 +sQz ∂2Qy(r,t)
∂z2 − ∂fL
∂Qy(r,t) ,
∂P(r,t)
∂t =L
sPx∂2P(r,t)
∂x2 +sPy∂2P(r,t)
∂y2 +szP∂2P(r,t)
∂z2 −∂(fL−fE)
∂P(r,t) . (8)
Again in practical calculations, we take the values of the parameters given by Artyukhinet al.unless stated otherwise,
as shown in Table I (including L=1) obtained from the first-principles calculations on YMnO3 [25]. The periodic
FIG. 2. The plotted ground-state solutions of parametersQ(a), P (b), andF(c) as a function ofErespectively, as calculated from Eq. (6). Here the gradient energyFGis not included. The dashed lines are the unstable solutions.
boundary conditions are employed and the lattice size used will be specifically stated below. The finite difference method is employed and the initial lattice is set as
Qx(r)=0, Qy(r)=0, P(r)=0. (9) Given the above setting, the domain structures under differentEwill be calculated upon sufficiently big number of iterations of Eq. (7). Here, an issue deserved for clarification is the time-dependent kinetic factor. This issue can be clearly revealed by tracking the ferroelectric hysteresis obtained at dif- ferent iteration steps. To go ahead, we setEas a sine function
FIG. 3. The calculated ferroelectric (P-E) hysteresis loops at different frequenciesωof the ac electric fieldE. The values ofω are inserted numerically.
of timet, i.e.,E=Emsin(2π ω·t) whereEm=0.8 andωis the frequency, and apply this field to a well-evolved topological structure, whose size is 2048x ×2048y ×1z with x=y=z=0.2 nm. TheP-E loops with differentω values as labeled are shown in Fig. 3. While remarkable dependence of the loop shape and area onω is presented, it is suggested that a frequencyω <10−5(1/iteration steps) is sufficiently low to guarantee a quasistatic dynamic evolution.
Hereafter, all the data to be presented are obtained from the simulations under a given constantEafter 105iteration steps (ISs), unless stated elsewhere.
III. RESULTS AND DISCUSSION
A. Evolution of domain structure and domain-wall motion In this simulation, a well-evolved topological structure in a cubic lattice of 512x × 512y × 512z with x = y=z=0.2 nm is used as the initial structure. For a better illustration, we cut a part of this lattice as shown in Fig. 4(a), and a set of successive domain patterns ob- tained atE=0.16,0.32, and 0.66, respectively, are shown in Figs.4(b)–4(d). Afterwards, the electric field is reset asE=0, and the eventual domain pattern is shown in Fig. 4(e). A comparison of Fig.4(a)with Fig.4(e)suggests that the domain structure is path dependent, a character of ferroelectricity, noting that the topological invariances of the two patterns remain the same.
TheE-dependent structure evolution exhibits several fea- tures. First, as E >0, those P↑ domains expand and those P↓ domains shrink. In the high-field range, one sees that those P↓ domains shrink into narrow stripes, more or less TABLE I. The values of physical parameters used in the present calculations (e.g., YMnO3) [23]. Note: These values are used in all the calculations unless stated elsewhere.
a0(eV ˚A−2) b0(eV ˚A−4) c0(eV ˚A−6) c0(eV ˚A−6) aP(eV ˚A−2) g(eV ˚A−4) g(eV ˚A−4)
−2.626 3.375 0.117 0.108 0.866 1.945 9.931
sxQ(eV) syQ(eV) szQ(eV) sxP(eV) syP(eV) szP(eV) L
5.14 5.14 15.4 8.88 8.88 52.7 1
FIG. 4. The simulated three-dimensional domain patterns at differentE. The magnitude ofP is scaled by the color bar. The varying sequence of electric fieldEis (a) 0.00→(b) 0.16→(c) 0.32→(d) 0.66→(e) 0.00. A small rectangle area covering a domain wall is marked for reference in Fig.5.
like wide “domain walls.” Second, as a signature of topo- logical protection, those domain walls may move but those vortex/antivortex cores do not [17,38,39]. Third, as mentioned above, the domain structure does not recover back to the initial state after the poling process. For a quantitative discussion, the domain-wall motion can be described by the spatial profiles of order parameters and energy density. These profiles taken from a small rectangle area as marked in Figs. 4(a)–4(c) are plotted in Figs. 5(a)–5(c). The line profiles of P, Q, and, and total-energy density F along the dashed lines in Figs.5(a)–5(c)are plotted in Figs.5(d)–5(g). The domain-wall width is roughly 0.5−1.0 nm, the magnitudes of bothP and Q in theP↑ domains slightly increase and those in theP↓ domains slightly decrease either with increasingE, the energy density F in the P↑ domains falls down but that in the P↓ domains rises up with increasingE, and parameterin the P↑ and P↓ domains remains unaffected by field E. These evolution features mark the transition from the type-I pattern to the type-II one, driven by electric field.
For the kinetics of domain-wall motion, we track the walls under different fields in a lattice with much bigger in-plane size (2048x × 2048y × 1z) and the data exhibit much better accuracy (x=y=z=0.02 nm). The perfect linear dependence of wall speedvonEis evaluated, as shown in Fig. 5(h), and can be confirmed by checking the energy difference (F↑−F↓) between the P↑ and P↓ domains as a function ofE, as shown in Fig.5(i). Clearly, the perfect linear dependence of (F↑−F↓) onE reasonably yields the linear v(E). It should be mentioned that this linear dependence is evaluated by ignoring the role of elastic energy which
may be non-negligible for some ferroelectrics, and therefore experimentally measured dependence for those systems may deviate from this linear relation.
It is noted that the above simulations are carried out in the low-field range. In this case, these domain walls move but none of them is destroyed by wall fusion. Here, we use the terminology “fusion” in the sense of the spatial resolution of ∼0.2 nm or so, to be discussed later. This low-field topological protection property, as observed experimentally too, has been exemplified to explain why the vortex-antivortex pair annihilation cannot occur. The microscopic mechanism is that the commensurability of partial unit-cell shifts across the paired walls must be satisfied, ensuring the topological protection and reflecting the topological invariance of the present structure [38]. Beyond this property, one may be more interested in checking the domain structure in the high-field range, which may not be realizable experimentally at the current stage due to the extremely high-field threshold, but can be virtually accessed in our simulations. In this case, if two neighboring walls merge, the two domains aside this newly born wall would have the same polarization orientation, suggesting the possibility of topological protection disabling.
Certainly, this emergent phenomenon, coined as the fusion of domain walls, is of special interest and will be discussed in this work.
B. Fusion of domain walls
We focus on a local area where two AP+FE domain walls with three neighboring domains are included, as shown in
FIG. 5. The stationaryPcontours of the rectangle area as marked in Fig.4, as obtained upon implication of different electric fieldE= 0.00 (a), 0.16 (b), and 0.32 (c). The spatial line profiles of parameters P(x),Q(x),(x), andF(x) across the domain wall, as marked in (a)–(c) are presented in (d)–(g) respectively. The domain-wall motion speedv(h) and energy density difference (F↑−F↓) (i) as a function of electric fieldErespectively.
Fig.6(a). The trimerization phase difference across the two walls here is 2=2π/3. As shown in Figs.6(a)–6(d)as a consecutive sequence, the enhanced electric field gradually enforces the in-between P↓ domain to shrink. Its two side walls come to meet and merge into one wall in the sufficiently high E=0.66. This wall also has the upward polarization which is smaller than that of the two side domains, as shown in Fig. 6(f) where the P profile along the dashed white line in Figs.6(a)–6(d)is plotted. This process is called the fusion of ferroelectric domain walls, corresponding to a transition from =π/3 to=2π/3 with increasing E. This fusion behavior is quite unusual which never occurs for normal ferroelectrics where such a field-driven domain merging is spontaneous. The underlying physics is naturally associated with the interlocking of the ferroelectric domains with antiphase structure, while this antiphase structure is robust against electric field unless the field is sufficiently high.
This fusion process is also illustrated by the spatial profiles of Q and as well as F, as plotted in Figs. 6(g)–6(i). It is clearly seen from Fig. 6(g)that the two valleys of Q(x), well separated atE=0, come to meet each other and merge into one valley above E∼0.32. Similar behavior can be seen from Fig. 6(i) for F(x) that has two peaks instead of
FIG. 6. The stationaryP contours across two ferroelectric do- main walls, as obtained upon sequential implication of different electric fieldsE=0.00 (a), 0.16 (b), 0.32 (c), 0.66 (d), and 0.00 (e). The spatial line profiles of parametersP(x),Q(x),(x), and F(x) across the domain walls, as marked in (a)–(e) are presented in (f)–(i) respectively.
valleys. The plateau between the two valleys (peaks) gradually disappears. In particular, theFpeak is quite large atE∼0.66, suggesting that this merged wall is highly unstable once the field is removed when this merged wall will again split into two separate walls, as shown by the profiles in Fig.6(e)whenEis back toE=0 in sequence. For profile(x), the fusion of the two walls can also be identified, corresponding to the transition from=π/3 to=2π/3. It should be mentioned that such a fusion has been neither observed experimentally nor predicted from the topological structure.
Besides the phenomenological analysis, it is helpful to look at the domain-wall fusion process from the viewpoint of atomic configuration, based on our simulated results. We consider two domain walls, each of which has a trimerization phase difference =π/3. It is assumed that a possible fusion of the two walls produces a relatively wide “wall”
with=2π/3, as described above. The two walls viewed on the ab plane, marked by two black dashed lines, are drawn in Fig.7. For convenience of discussion, the two walls in the fusion state are assumed to be close to each other, constituting the wide wall with=2π/3, in other words, the middle wall with downward polarization is sandwiched by two P↑ domains. Based on existing experimental and theoretical results, the wide wall with =2π/3 can be constructed via two different manners [17,40] shown in Figs. 7(a) and 7(b) respectively, taking hexagonal YMnO3 as an example. In Fig.7(a), the neighboring Y ions take the (half-up)-down-down-(half-up) alignment across the “wall”
while in Fig.7(b)the Y ions take the down-up-down alignment.
In the absence of electric field, the polarization inside the
“wall” is downward (marked by the violet rectangles) due to
FIG. 7. The lattice configurations for two types of AP-only domain walls with=2π/3. The atomic structures of type-A (a) and type-B (b) domain walls seen along the [001] and [100] axes. The “up” and “down” indicate the directions of the Y ions’ displacements from the paraelectricP63/mmcto the ferroelectricP63cmstructure. The arrow in the in-plane lattice configurations indicates the trimerization phase orientations. The coarse open arrow indicates the polarization which changes its sign from negative value to positive value with increasing electric field.
the down-shifted Y ions. The electric field along thec axis has twofold impacts. On one hand, the increasing electric field would push the Y ion layers upward, leading to the variation of the Y-O bond length. The polarization inside the “wall”
would change from negative value to positive value gradually.
On the other hand, the separation between the down-shifted and up-shifted Y-ions inside the “wall” becomes smaller with larger field, as marked by the blue rectangles. This reduced separation certainly weakens the trimerization inside the wall.
Our simulations show that the Q value in the wall with =2π/3 does decrease with increasingE, implying that the inclination magnitude of MO5triangular bipyramids inside the wall becomes small. Both scenarios shown in Figs. 7(a) and7(b)are consistent with our simulations shown in Fig.6.
Here it should be mentioned that such a fusion of two neighboring AP+FE domain walls seems to be a transient state driven by a sufficiently high electric field, which is certainly unstable at E=0. The fused domain wall would decompose into two walls when the electric field returns back toE=0, as illustrated in Figs.6(d)and6(e). Indeed, neither such a fusion of two neighboring walls with=π/3 into one AP-only wall with=2π/3 nor such a decomposition of one AP-only wall with=2π/3 into two AP+FE walls with=π/3 have been observed experimentally. It is thus misleadingly believed that a domain wall with=π/3 is so robust that a fusion of two such domain walls into one wall with=2π/3 is impossible. Nevertheless, so far all the experimental observations were performed under E=0 after the high-field poling treatments, and the wide wall with =2π/3, which can be produced in the high-field range, becomes improbably observed atE=0. Our simulations just reveal that such a wide wall is indeed generated by the high electric field. Furthermore, such fusion and decomposition
events do not break the topological invariance of the domain structure, and for all cases the variation of parameter counting around any vortex or antivortex core remains to be 2π, no matter whether the walls are with either=π/3 or =2π/3.
C. Evolution of vortex-antivortex pairs
The simulations presented in the above sections deal with the domain walls. Now we discuss the evolution of vor- tex/antivortex cores in response to increasing electric field. For unveiling the general characters of the cores, our simulations are performed in a grid of 2048x × 2048y × 1zwith x=y=z=0.02 nm, noting that the gridding is much finer than the previous cases so that every particular of the core region can be tracked to avoid any skeptical conclusion.
A set of electric fields along thecaxis are applied step by step, and 105iterations are cycled for reaching the quasistatic domain structure at each step. For convenience of discussion, we present the results on one vortex core without losing the generality. The typical spatial profiles of parametersP,Q,, andF, at different fields are shown in Fig.8forE=0,0.20, and 0.66, respectively.
The core region atE=0 is mapped in Fig.8(a)by these parameters. First, the P contour shows the normal vortex pattern. Second, a conicalQpattern in this region is obtained.
Outside this region, one sees the hexapetalous Q pattern, mapping the six interlocked AP-FE domains. It is found that Q=0 and P ∼0 at the core tip, implying that the meeting point of the six AP-FE domains accommodates a high symmetric phase (P ∼0 and Q=0), consistent with the Kibble mechanism that the defects are remnants of the parent phase trapped within the lower symmetry phase
FIG. 8. The spatial contours of parametersP,Q,, andFover a vortex core region, given different electric fields (a)E=0, (b) E=0.20, and (c)E=0.66, respectively.
[41]. The six domains are separated by a phase difference =π/3 and the core tip exhibits the highest energyF. AsE >0, e.g., atE=0.20, as shown in Fig.8(b), the three P↓ domains shrink in size (width) in compensation with the three expandingP↑ domains. TheQpattern around the core tip remains nearly unchanged, but the contour outside this region becomes trivalvelike, resulting from the fusion of two AP+FE domain walls with=π/3. This fusion can be also clearly seen from thepattern andF pattern, andF at the core tip is the largest. Obviously, the incredible robustness of the vortex/antivortex cores against the electric field is the consequence of the topological protection.
AsEfurther increases up to 0.66, as shown in Fig.8(c), the parameter contours over the core region change a lot. TheP at the core tip, which should be very small atE=0, increases gradually withEin a linear manner (Fig.9). The core region, as circled in Fig.8(c1), becomes ferroelectric (P =0). In fact, the whole lattice hasP >0 everywhere and is occupied by a singleP↑ domain embedded with some stripelike regions of smallerP↑. ThisE∼0.66 is the critical field beyond which the sixfold ferroelectric domain structure associated with the vortex/antivortex cores is destructed. Although the sixfold AP+FE domain structure is transformed into a threefold AP-only domain structure duo to the fusion of two AP+FE domain walls into an AP-only domain wall, the topological property of structural antiphase domain patterns and the Q contour over the core region remain unchanged atE=0.66.
In fact, what happens at the vortex cores is the decoupling between FE domains and AP domains, comparing the P, Q, and patterns under E=0 and E=0.66. Another
FIG. 9. The calculated polarization P at the vortex core tip, energy densityFat the vortex core tip, and difference (F) in energy density between the vortex core tip and outer domain, as a function of electric fieldErespectively.
interesting fact is that the energy densityF around the core tip gradually decreases with increasing E, indicating that the vortex structure evolution is spontaneous. In particular, the energy difference F between the vortex core tip and the domains stemming from this tip also decreases with increasingE, as shown in Fig.9.
The above outlined results indicate that the sixfold FE domain structure can be destructed by electric field, and this further raises a question: Will the topological invariance prop- erty associated with the structural vortex-antivortex domain patterns be broken by sufficiently high field? However, here the topological invariance in terms of the trimerization phase around the vortex core remains to be 2π, upon a field up to E∼0.66, i.e., this invariance property is maintained. Next, it is interesting to check this property by further increasingE.
D. Electric field induced phase transition
We check the stability of this threefold antiphase structure associated with a vortex upon further increasingE. We first check theQ=0 property at the core tip and our simulations show that this property can be maintained untilE=1.48 from 0.66. Thein situtracking of the spatial contours of parameters P, Q, and F at E=1.48 are carried out and the snapshot patterns of them after cycling for 6×105ISs are plotted in Fig.10(a). It is shown that the whole core region has nonzeroQ values except the core tip point at whichQ=0. TheP contour over this region is positive and its values around the core tip are larger than those outside. The energy density F along the antiphase walls is higher than that inside the antiphase domains.
It is surprising to observe thatE=1.48 is a critical electric field. When the sixfold ferroelectric domain structure over the vortex core is already destructed at a field ∼0.66, this critical fieldE=1.48 marks the destruction of the threefold antiphase structure over the core region. If one increases slightly the field toE=1.49 from 1.48, a substantial change of the P, Q, and F spatial contours is identified. For a better illustration, we show in Fig.10(b) the contours after the cycling for 5×104ISs atE=1.49. The contours after
FIG. 10. The spatial contours of parametersP,Q, andF over a vortex core region, given different electric fields (a)E=1.48 after cycling for 6×105ISs, (b)E=1.49 after cycling for 5×104ISs, and (c)E=1.49 after cycling for 8×104ISs, respectively.
the cycling for 8×104ISs are presented in Fig.10(c)for a comparison. Several features deserve for highlighting here.
First, the Qcontour over the whole region as E <1.49 is nonzero except at the core tip, implying that the topological invariance property is still maintained. However, atE=1.49, thisQ=0 point expands rapidly with time, as shown by the snapshot contours in Figs. 10(b2)and10(c2)at 5×104ISs and 8×104ISs respectively. The appearance of such aQ=0 region rather than an isolated point at the tip is a clear mark of the topological invariance missing. In this sense, it is meaningless even that the variance ofsurrounding the core tip is 2π. In other words, a lattice withQ=0 no longer has the trimerization property and thus the topological invariance.
Second, polarization P over this region increases rapidly in thisQ=0 region, as shown in Figs.10(b1)and10(c1). The magnitude ofP in this region is much larger than that outside.
It is found in our calculations that the dependence ofP onE in theQ=0 region still follows the relationshipP =E/aP, whereaPis the proportional coefficient. Third, this process is energetically favored, as shown by the much lower energy F in this region at E=1.49 than that at E=1.48. This behavior seems to be quite consistent with the scenario of a field-induced phase transition, marking the disappearance of the topological invariance. The details of this transition will be discussed elsewhere in the future.
This unusual missing of the topological invariance certainly needs its lattice configuration consistence, and here we discuss the possible evolution of lattice structure corresponding to this phase transition. It is understood that there are three structural antiphase domains (α+,β+, andγ+) around a vor- tex/antivortex core. A strong competition among theα+,β+, andγ+domains upon sufficiently high fieldEmay lead to the domain decoupling at the core. In addition, it is revealed by
FIG. 11. The side view of the paraelectric (PE,P63/mmc) (a) and ferroelectric (FE,P63cm) (b) crystal structures of hexagonal YMnO3. The “up” and “down” indicate the directions of the Y ions’s displacements from the PE to the FE structure. Distinct Y-Oepbonds resulting from the inhomogeneous Y ions distortion are highlighted with green and red lines, respectively.
Fig.2(a)that theQvalue inside theα+,β+, andγ+domains only varies slightly with increasingE. This implies that the electric field induced phase transition occurring at the core is not induced by the competition among theα+,β+, and γ+ domains.
A possible microscopic origin for this phase transition can be proposed. For the paraelectricP63/mmc phase, the lengths of two neighboring Y-Oepbonds along thecaxis are equivalent, as shown in Fig. 11(a). The shifting of Y ions along thecaxis during the structural phase transition from the centrosymmetricP63/mmcto the ferroelectricP63cmmakes the two Y-Oepbonds nonequivalent: one becomes longer and the other shorter, leading to a net electric polarization. Figure 11(b) shows the ferroelectric P63cm structure, where the polarization direction is along thecaxis. A positive electric field favors the red Y-Oep-Y bonding rather than the green Y-Oep-Y bonding. This is probably the reason why parameter Qin theP↑ states varies slightly with increasingE >0 but that in theP↓states varies remarkably, as shown in Fig.2(a).
Therefore, the positive electric field more strongly favors the crystal structure in which all Y-Oep-Y bonding aligns as the red bonds indicate. In this case, all the Y ions on the same layer no longer show intershifting from each other along thecaxis.
Therefore, no more tilting of the MnO5triangular bipyramids occurs, i.e., no more trimerization (Q=0). This suggests that the trimerization-induced ferroelectric phase is stable only in the low-field range, and it will be replaced with another ferroelectric phase withQ=0 in the sufficiently high electric field range where the topological invariance is no longer a valid property. This electric field induced transition from the trimerization-induced ferroelectric phase to an electric field induced ferroelectric phase with Q=0 occurs at the vortex/antivortex core.
E. Discussion
To this end, we have focused on the domain-wall fusion and vortex/antivortex core destabilization as a signature of the topological invariance disabling. The results are indeed interesting and informative for understanding the underlying physics. However, we have shown how the topological domain
structure evolves in response to electric field along thecaxis, but no case of electric field imposed along the other direction is discussed. It is believed that a low in-plane electric field would impose no effect on the ferroelectric domain structure, the associated piezoelectric and electrostrictive effects, as well as electric static effect may distort the lattice structure, weakening or strengthening the lattice trimerization. This consequence if any cannot be discussed in the present framework of phenomenological Ginzburg-Landau theory. However, thec- axis electric field as a case without losing the generality does reveal the major features of domain structure evolution.
On the other hand, a major drawback of this understanding is that all the conclusions are only theoretical predictions based on the proposed Ginzburg-Landau theory on hexagonal manganites RMnO3, and no evidence with the applicability of this theory in the high electric field has been available.
This uncertainty can be partially cancelled by our proposed microscopic models for the main results, as shown in Figs.7 and11. Second, no experimental observation on any of the main conclusions in this work is available. This work thus appeals for experimental check, noting that the TEM and other domain structure probing techniques available so far can be improved to operate under certain electric field applied to the samples.
Furthermore, since extremely large electric fields are applied in this work, higher-order terms in the Landau free-energy expansion may become non-negligible and the consequent result may be different. In order to check this point, the fourth-order term inP(P4), is added to the Landau energy FLof Eq. (2). Now, the Landau energy has the following form:
FL = a0
2Q2+b0
4Q4+c0
6Q6+c0
6Q6cos 6
−gQ3Pcos 3+g
2Q2P2+aP
2 P2+bP 4 P4,
(10) wherebP is the precoefficient to the fourth-order term inP. Based on the standard procedure, Eqs. (5) and (6) now can be written into Eqs. (11) and (12) shown below, respectively:
a0Q+b0Q3+(c0+c0)Q5−3gQ2Pcos 3 +gQP2=0, −gQ3cos 3+gQ2P +aPP +bPP3=0, =0,±π/3,±2π/3,π, (11) a0Q+b0Q3+(c0+c0)Q5−3gQ2P +gQP2 =0,
−gQ3+gQ2P +aPP +bPP3 =0, =0,±2π/3, P↑,P >0, a0Q+b0Q3+(c0+c0)Q5+3gQ2P +gQP2 =0, gQ3+gQ2P +aPP +bPP3 =0,
= ±π/3,π, P↓,P <0, (12)
where all coefficients exceptbPin Eq. (10) have been obtained from the first-principles calculations on YMnO3and are listed in TableI.
FIG. 12. The plotted ground-state solutions of parametersQ(a), P(b), andF(c) as a function ofErespectively when the fourth-order term inP is added to the Landau energy. Along the violet arrow, the curves successively correspond to the results atbP =0.0,0.8,1.6,2.4, and 10. The red arrow marks the shifting of the intersection point of the twoF(E) curves for the P↑ state and theQ=0 state as bP
increases.
The value of coefficientbP must be available for a new calculation. Obtaining this value from the first-principles calculations is beyond the scope of this work. Instead, we take a set ofbP values (bP =0,0.8,1.6,2.4, and 10) to perform the calculation and see what is the effect of its inclusion. The results are plotted in Fig.12, where the fine and long arrows in each plot indicate increasedbP. It should be mentioned that bP =0 corresponds to the case of no inclusion of the fourth-order term inP.
Several features can be highlighted regarding the inclusion of the fourth-order term in P. First, an inclusion of this term does not change the qualitative behaviors of Q(E), P(E), and F(E) for both theP↑ state and theP↓ state, as shown in Figs.12(a)–12(c)of the paper, while the quantitative dependence is yet insufficient. No big qualitative difference can be identified even whenbP is very large, indicating that the role of the fourth-order term in P is very important. A remarkable difference appears only in extremely large E.
Second, for the state ofQ=0, the rapid decreasing tendency of free energy F with increasing E is weakened when bP
increases, as shown in Fig.12(f). Third, the intersection point of the twoF(E) curves for theP↑state and theQ=0 state, where the two states have equal energy, moves towards a large Evalue whenbPincreases, as shown in Fig.12(f), and so does the intersection point of the twoF(E) curves for theP↓state and theQ=0 state. This behavior implies that a biggerbPwill
induce a widerEwindow for the trimerization state in which the trimerization state has the lowest free energy and higher stability. A destabilization of the trimerization state requires a largerE.
Another point is the huge gap in magnitude of electric field between our simulations and experiments. While this gap implies the necessisity to improve the Ginzburg-Landau theory, the main results as simulated in this work would make sense qualitatively. In fact, the electric field needed in practical experiments can be roughly estimated based on our results.
For instance, the electric field needed to prompt the fusion of two AP+FE domain walls into one AP-only wall for measurements can be∼107 kV/cm since the coercive field in this work is∼0.247, and the electric field being able to disable the topological antiphase domain structure is∼241 kV/cm.
In this sense, this prediction can be confirmed experimentally.
IV. CONCLUSIONS
In this work, we have investigated in a systematic way the responses and relevant dynamics of the real-space topological domain structure in hexagonal manganites RMnO3 (e.g., YMnO3) against external electric field along thecaxis, by us- ing the phase-field simulations based on the phenomenological Ginzburg-Landau theory. Our simulations have led to several major consequences. First, it has been revealed that in the low-field range the topological structure is robust against the
electric field, characterized by the immobile vortex-antivortex core while the ferroelectric domain pattern is driven from the type-I pattern to the type-II one. Second, it has been identified that a high electric field can make the two interlocked AP+FE domain walls with=π/3 fuse into one AP-only wall with =2π/3, by the field-driven disappearing of theP↓ domain sandwiched by the two interlocked AP+FE walls. On the other hand, the high electric field can also seriously deform the topological vortex/antivortex structure, characterized by the disabled sixfold ferroelectric domain pattern and the reserved topological invariance in terms of the structural trimerization. This topological invariance can be further destructed by sufficiently high electric field which sup- presses the structural trimerization, and the vortex/antivortex core region can be replaced by a trimerization-free state (Q= 0), which seems to be a phase-transition-like process. The present simulations have revealed a set of topological domain structures in response to external electric field, representing a substantial forward step in our understanding of real-space topological structure in hexagonal manganites.
ACKNOWLEDGMENTS
This work was supported by the National Key Research Program of China (Grant No. 2016YFA0300101) and the National Natural Science Foundation of China (Grants No.
51431006, No. 11374147, No. 11234005, No. 11504048, and No. 11774106).
[1] N. D. Mermin, The topological theory of defects in ordered media,Rev. Mod. Phys.51,591(1979).
[2] W. H. Zurek, Cosmological experiments in superfluid helium?
Nature (London)317,505(1985).
[3] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Topological semimetal and Fermi-arc surface states in the electronic structure of pyrochlore iridates,Phys. Rev. B 83, 205101(2011).
[4] Y. L. Chen, J. G. Analytis, J. H. Chu, Z. K. Liu, S. K. Mo, X. L.
Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R.
Fisher, Z. Hussain, and Z. X. Shen, Experimental realization of a three-dimensional topological insulator, Bi2Te3,Science325, 178(2009).
[5] Y. Du, F. Tang, D. Wang, L. Sheng, E.-j. Kan, C.-G. Duan, S. Y.
Savrasov, and X. Wan, CaTe: A new topological node-line and Dirac semimetal,npj Quantum Mater.2,3(2017).
[6] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems,J. Phys. C6,1181 (1973).
[7] D. H. Van Winkle and N. A. Clark, Direct measurement of orientation correlations in a two-dimensional liquid-crystal system,Phys. Rev. A38,1573(1988).
[8] O. V. Yazyev and S. G. Louie, Topological defects in graphene:
dislocations and grain boundaries, Phys. Rev. B 81, 195420 (2010).
[9] K. Suenaga, H. Wakabayashi, M. Koshino, Y. Sato, K. Urita, and S. Iijima, Imaging active topological defects in carbon nanotubes,Nat. Nanotechnol.2,358(2007).
[10] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Real-space observation
of a two-dimensional skyrmion crystal,Nature (London)465, 901(2010).
[11] F.-T. Huang, B. Gao, J.-W. Kim, X. Luo, Y. Wang, M.-W.
Chu, C.-K. Chang, H.-S. Sheu, and S.-W. Cheong, Topological defects at octahedral tilting plethora in bi-layered perovskites, npj Quantum Mater.1,16017(2016).
[12] Y. Li, Y. Jin, X. Lu, J.-C. Yang, Y.-H. Chu, F. Huang, J. Zhu, and S.-W. Cheong, Rewritable ferroelectric vortex pairs in BiFeO3, npj Quantum Mater.2,43(2017).
[13] Y. Geng, H. Das, A. L. Wysocki, X. Wang, S. W. Cheong, M.
Mostovoy, C. J. Fennie, and W. Wu, Direct visualization of magnetoelectric domains,Nat. Mater.13,163(2014).
[14] H. Das, A. L. Wysocki, Y. Geng, W. Wu, and C. J. Fennie, Bulk magnetoelectricity in the hexagonal manganites and ferrites, Nat. Commun.5,2998(2014).
[15] S. Dong, J.-M. Liu, S.-W. Cheong, and Z. Ren, Multiferroic ma- terials and magnetoelectric physics: Symmetry, entanglement, excitation, and topology,Adv. Phys.64,519(2015).
[16] S. C. Chae, Y. Horibe, D. Y. Jeong, S. Rodan, N. Lee, and S.-W.
Cheong, Self-organization, condensation, and annihilation of topological vortices and antivortices in a multiferroic,Proc. Natl.
Acad. Sci. U.S.A.107,21366(2010).
[17] T. Choi, Y. Horibe, H. T. Yi, Y. J. Choi, W. Wu, and S.-W. Cheong, Insulating interlocked ferroelectric and structural antiphase domain walls in multiferroic YMnO3,Nat. Mater.9, 253(2010).
[18] C. J. Fennie and K. M. Rabe, Ferroelectric transition in YMnO3 from first principles,Phys. Rev. B72,100103(2005).
[19] T. Lonkai, D. G. Tomuta, U. Amann, J. Ihringer, R. W. A.
Hendrikx, D. M. Többens, and J. A. Mydosh, Development of
