Homotopy-Theoretic Study & Atomic-Scale Observation of Vortex Domains in Hexagonal Manganites

Essential structural properties of the non-trivial “string-wall-bounded” topological defects in hexagonal manganites are studied through homotopy group theory and spherical aberration-corrected scanning transmission electron microscopy. The appearance of a “string-wall-bounded” configuration in RMnO3 is shown to be strongly linked with the transformation of the degeneracy space. The defect core regions (~50 Å) mainly adopt the continuous U(1) symmetry of the high-temperature phase, which is essential for the formation and proliferation of vortices. Direct visualization of vortex strings at atomic scale provides insight into the mechanisms and macro-behavior of topological defects in crystalline materials.

Movie II | Simulated annealing process with T a = 1250 K. Movie IIamplitude.mov and Movie IIphase.mov show the distribution of order parameter field Q and φ, respectively. The temperature is high enough to drive system into U (1) symmetry. This provides the necessary condition for vortex core formation.
Movie III | Simulated annealing process with T a = 1300 K. Movie IIIamplitude.mov and Movie IIIphase.mov show the distribution of order parameter field Q and φ, respectively. The temperature is high enough to drive system into U (1) symmetry. This provides the necessary condition for vortex core formation.

Supplemental Information
Section S1. Landau free energy and degeneracy order parameter space of RMnO 3 For hexagonal manganites RMnO 3 , S. Artyukhin et al. have given the free energy density expansion based on Landau theory [1]. It can be represented in the following form: where f bulk is the bulk free energy density and f grad is the gradient energy density. Since the primary order parameter of improper ferroelectric transition has at least two components, scalar variables Q and φ are used to describe the amplitude and azimuthal angle of trimerizing tilt, respectively. Spontaneous polarization is a secondary order parameter in this system, which is quantified by the amplitude (P) of polar Γ 2 mode. It is triggered by the distortion of MnO 5 3+ bipyramids. All coefficients in equation (1) had been determined by a series of first-principles calculations in Ref. 1. The coefficient of Q 2 varies with temperature T, T s ≈ 1270K is the critical temperature of structural phase transition for YMnO 3 . The bulk free energy density can be decomposed into two parts (separated by blankets in equation (1)), denoted as f 1 and f 2 , respectively, where f 2 breaks the full symmetry of f 1 .
Considering a crystal in mono-domain state, the equilibr ium state of order parameters can be determined by minimization of f bulk with respect to Q, φ and P. This gives Q 0 = 0.956 Å (when T = 0 K), φ belongs to Φ = [0, π/3, 2π/3, π, 4π/3 and 5π/3], and = 3 cos 3φ/( P + ′ ). Hence, the order parameter space R 1 of this uniform system is six evenly spaced points on the edge of a circle. As the temperature rises, both the bulk free energy density and relative weight of f 1 in f bulk increases, as shown is Fig. S1. When temperature approaches T s , f bulk varies little with different φ and f 1 becomes dominant in f bulk (i.e. |f 1 | >> |f 2 |), the minimization of f bulk gives decreased Q value and depends little on the variation of φ. Then the degeneracy space expands from R 1 to R 2 .
When T > T s , -a(T -T s )/2T s becomes positive, and Q = 0 gives the minimum of f bulk . Then the circle R 2 shrinks into a single point R 3 . Above all, the degeneracy order parameter space is a function of temperature.

Section S2. Homotopy groups and exact homotopy sequences
According to homotopy theory, topological defects can be classified by the elements of homotopy groups associated with the symmetry of order parameter space. The zeroth homotopy group π 0 (R, x 0 ) classifies two dimensional domain walls (x 0 is a base point in order parameter space R), the fundamental homotopy group π 1 (R, x 0 ) classifies one dimensional vortex lines, and the second homotopy group π 2 (R, x 0 ) classifies the point defects. In these groups, x 0 is a base point in order parameter space R and it can be omitted when R is a connected space. These groups are the special cases of the so-called relative homotopy groups. If A is a subset of R which contains the base point x 0 , the nth relative homotopy group can be expressed as π n (R, A, x 0 ). When the degeneracy parameter is restricted to A, for instance on surfaces, or in regions of small parameter variations, we can use relative homotopy groups to classify these topological defect structures. In addition, relative homotopy groups have the following relationship: … → π +1 ( , , 0 ) → π ( , 0 ) → π ( , 0 ) → π ( , , 0 ) … …π 2 ( , , 0 ) → π 1 ( , 0 ) → π 1 ( , 0 ) → π 1 ( , , 0 ) → π 0 ( ) → π 0 ( ), in which "→" represents homomorphic mapping, and this series of mappings is called exact homotopy sequence. In this sequence, the image (im) of one homomorphism always equals to the kernel (ker) of the next homomorphism, i.e. im (i) = ker (j). after an annealing process above T s , the coexistence state of vortices ("string-wall" bounded type) and circle domains can be observed.

Section S3. Evolution of order parameter field with temperature and external electric field
To consolidate the result obtained from homotopy theory analysis, the evolution of order parameter field with temperature is simulated on a two dimensional rhombic lattice with periodic boundary condition, based on the Landau phenomenological model (see equation 1). Each site of the lattice is assigned two variables Q and φ whose values are updated by the Metropolis algorithm. Three typical snapshots of the evolution are shown in Figs. 2(a-c), corresponding to three states at T > T s , T s > T >> 0 and T = 0, respectively. Figure 2(a) shows a defect-free homogenous state with all Q values in close proximity to 0. Just below the phase transition temperature, Q increases slightly in the demand of minimization of bulk free energy which is still non-sensitive to the value of φ, so there is no preference for φ. However, the spatially configurations of order parameter field has noticeable impact on the gradient energy, so smooth variation of φ among adjacent sites is energetically preferred. The system in this state is in analogy with the low -temperature phase of the x-y model in spin system because their degeneracy spaces both adopt U(1) symmetry, then topological excitations (bounded vortexantivortex pairs without domain walls) will appear to minimize the total energy, as shown in Fig. 2 In order to reduce gradient energy, some new order parameters which beyond the degenerated space may appear in spite of the additional cost in bulk free energy. Close to the wall, the increase of gradient energy density results in deviation of Q and φ from R 1 . The value of φ is no longer confined to those six distinct values as in the region away from defects. Instead, it spreads into continuous distribution. Supposing that φ changes sharply between two adjacent domains, such as from 0 to π/3 between  + and  -, then the gradient energy diverges. In order to remove this singular gradient, φ has to vary from 0 to π/3 continuously in the order parameter field even though this area will be narrow. This is coincided with the first principle calculation Having obtained the distribution of Q and φ, it is convenient to calculate the configuration of secondary order parameter P near vortex core by using equation = 3 cos 3φ/( P + ′ ). The result is shown in Fig.   S3(a), P changes sign in adjacent domains and varies to nearly zero inside vortex cores. Since P is proportional to the spontaneous polarization, the distribution of polarization can also be reflected from this figure. The polarization within a vortex core is obviously weaker than domain regions, so the stability of core should be less sensitive to external electric field. We also performed numeric simulations by introducing a static electric

Section S4. Relationship between order parameters and structural distortion
In the above analysis, we consider the RMnO 3 crystal as a continuous medium which is characterized by the order parameter field. In practice, the variation of order parameters can reflects the atomic distortion of crystal lattice. According to group theory, the condensation of the K 3 mode lowers the symmetry from P6 3 /mmc to P6 3 cm. Q reflects the relative strength of K 3 mode and its value is determined by the atomic displacement driven by condensation of K 3 mode, i.e. = √2( 3 Y1 ) 2 + 4( 3 Y1 ) 2 + 6( 3 Mn ) 2 + 6( 3 O1 ) 2 + 6( 3 O2 ) 2 + 2( 3 O3 ) 2 + 4( 3 O4 ) 2 Normally, the displacement value at each site decreases with smaller Q, so the distortion in R layers and tilting of MnO 5 3+ bipyramids changes accordingly with Q. Group theoretical analysis shows that there are three isotropy subgroups P6 3 cm, 3 ̅ 1 and P3c1 for the K 3 irreducible representation of space group P6 3 /mmc, so these three kinds of structures are symmetry-permitted to exist in RMnO 3 system, see Table S1. It is known that when φ∈Φ, the space group is P6 3 cm and R layers adopt the "up-up-down" or "down-down-up" configuration ( Fig. S5(a)). For space group 3 ̅ 1, the azimuthal angle φ of MnO 5 3+ bipyramids belongs to the set Φ' = {π/6, π/2, 5π/6, 7π/6, 3π/2 and 11π/6}. In such unit cell, one third of the R atoms distort upward, one third of them distort downward and the rest keep in position, as shown in Fig. S5(b). This configuration is φ is shown schematically in Fig. 3(b).
In order to verify this relationship, we carried out first-principles calculations on YMnO 3 using the projector-augmented wave method as implemented in Vienna Ab initio simulation package (VSAP). We used the LSDA+U method to describe the exchange-correlation functional with values of U = 8 eV and J = 0.88 eV that were the same as previous reports. The ground state structure was adopted from x-ray diffraction experiments and A-type magnetic configurations were adopted. All atoms were fully relaxed until the forces converged to 0.005 eV·Å . We constructed a structural model by rotating the tilting angle of the MnO 5 3+ bipyramids by π/6 and performed the structural relaxations with the positions of the apical oxygen atoms fixed and all other atoms relaxed. After the relaxation, we found that the configuration of Y atoms changes from "updown-down" to "up-middle-down". Though the middle Y atom is not accurately in the center of the nearest two Y atoms in the plane (may be induced by the fixed Mn-O ap bond length), a significant shift along the c axis occurs. As shown in Fig. S5(a) also indicate that φ deviates a little from Φ within several unit cells near domain walls. However, the impact of non-uniform order parameter field on the atomic structure should be more significant in the core regions because of higher gradient energy density. Continuous variation of φ is accompanied by remarkable decrease of Q around the cores, so the atomic structure of vortex core is more likely to be constructed by 3 ̅ 1 and P3c1, and the buckling state of R layers may present some differences from the normal P6 3 cm phase.

Section S5. Atomic structure of vortex core observed by STEM
Our experimental results have showed that the light substitution of In for Y could visibly increase the density of vortex cores, which is technically important for high-resolution TEM observations to image a vortex core within an area as small as a micron [3]. A typical atomic-resolved HAADF image of vortex core is shown is Fig. S6.

Section S6. Three-dimensional structure of ferroelectric domain in RMnO 3
The three-dimensional character of defects is an essential factor to consider in analyzing the twodimensional projected HAADF images. To visualize the distribution of vortex strings and domain walls, a 3D numerical simulation is performed on a 200×200×50 lattice with periodic boundary condition. The result obtained after a simulated annealing process is shown in Fig. S8(a). From the distribution of φ and Q, vortexantivortex pairs can be observed on any surface, this is consistent with the experimental observation as shown in Fig. S8(d). In order to visualize the distribution of defects alone, the stereostructure of a single domain and vortex strings are shown in Fig. S8(b) and (c), respectively. Different from the plate-like domain walls in BaTiO 3 or other proper ferroelectrics, the domain walls in RMnO 3 appear as curved surfaces, terminating at the surfaces or vortex strings. Two reasons may result in this difference: firstly, the annihilation of vortexantivortex pairs is an "astronomical" time scale process, so these curved walls are topologically protected by the vortices from high temperature; secondly, the YMnO 3 is a kind of semiconductor-like ferroelectric with relative high density of carriers, so the unfavorable electrostatic at the charged walls can be lowered by the directional movement of carriers. Similar to other string-type topological defects, such as disclinations, they form loops inside the medium or terminate on the surfaces. According to these features, we cannot expect that the walls are always parallel to the observation directions. Actually, the bend of domain wall can happen at the atomic level [2], and it is estimated that there are about 50 atoms in each column along the observation direction even within a 30 nm thick sample, so the overlap of domains is nearly inevitable.
Section S7. The "overlapping effect" near vortex core We construct a simplified three-dimensional atomic model which contains sharp domain walls and straight vortex string to analyze the impact of "overlapping effect" on the atomic contrast as shown in Fig. S9.

Section S8. Two-dimensional Gaussian fitting
Firstly, the image shown in Fig. 4(a) is divided into many conterminous rectangle pieces. In each piece, just one atom image is contained and the edges are defined by the minimum of local intens ity. Then, the intensity distribution in each piece is fitted by function: .
The intensity centers can be determined by the fitting parameter (x 0 , y 0 ) and the degree of ambiguous can be quantified by the longitudinal standard deviation σ y .

Section S9. Measurement of vortex core size
To provide a more intuitive result, we pick a region which cross the green box shown in Fig. 4(a), and a line scan is carried out for a single Y layer as shown in Fig. S11(a). The standard deviation and relative displacement along the c axis of Y 1 , Y 2 and Y 3 are plotted in Fig. S11(b) and (c), respectively. Since the standard deviation of atom Y 1 stays at a low level, no significant elongated contrast appear even across the core region. Meanwhile, the contrast at Y 2 site only becomes ambiguous in the vicinity of center and its contrast recovers with increasing distance; atom Y 3 shows ambiguous contrast all over within this region. The clear contrast at the Y 1 site shows that these atoms shift upward gradually from the lower site to a higher site by spanning over about 10 unit cells. Then the size of vortex core is estimated to be 50 Å. Fig. S1. Variation of f 1 /f bulk (solid lines) and f bulk (dotted lines) with temperature when φ is fixed to be 0 and π/6, respectively. With ris ing temperature, the ratio f 1 /f bulk keeps increasing except some oscillations near T s , this reveals that the magnitude of f 1 is always larger than f 2 , especially just below T s . The bulk free energy density f bulk also increases with temperature, and it becomes almost φ independent at high temperature region.   Tab. S1. Three isotropy subgroups for the K 3 irreducible representation of space group P6 3 /mmc.       Fig. 4(a). In each layer, we label the Y atoms as Y 1 , Y 2 , and Y 3 periodically from left to right. In order to analysis the atomic contrast in the red box labeled in (a), we fitted this Y layer by two-dimensional Gaussian function as shown in (b) for the variation of standard deviations and (c) for the shift of intensity centers.