Higher-order topological solitonic insulators

Pursuing topological phase and matter in a variety of systems is one central issue in current physical sciences and engineering. Motivated by the recent experimental observation of corner states in acoustic and photonic structures, we theoretically study the dipolar-coupled gyration motion of magnetic solitons on the two-dimensional breathing kagome lattice. We calculate the phase diagram and predict both the Tamm-Shockley edge modes and the second-order corner states when the ratio between alternate lattice constants is greater than a critical value. We show that the emerging corner states are topologically robust against both structure defects and moderate disorders. Micromagnetic simulations are implemented to verify the theoretical predictions with an excellent agreement. Our results pave the way for investigating higher-order topological insulators based on magnetic solitons.

Spin waves (or magnons) and magnetic solitons are two important excitations in magnetic system. Various topological states have been predicted in magnonic materials, such as topological magnon insulators [22][23][24][25][26] and magnonic Weyl semimetals [27][28][29]. Typical magnetic solitons include vortices [30,31], bubbles [32,33], and skyrmions [34][35][36][37], which are long-term topics in condensed matter physics for their interesting dynamics and promising applications [38,39]. It has been shown that the collective gyration motion of magnetic solitons exhibits the behavior of waves [40][41][42][43][44]. Recently, the periodic arrangement of ferromagnetic nanodisks, called magnonic crystal, has received much attention in the context of band structure engineering [43,[45][46][47]. When embracing the topological properties of vortex states, it paves the way for robust spintronic information processing. Topological chiral edge states are discovered in a two-dimensional honeycomb lattice of magnetic solitons [48,49]. However, all these topological magnonic and solitonic states are first-order in nature, according to the classification of topological insulators mentioned above. In this article, we predict a new class of higher-order topological insulators from the dynamics of magnetic solitons on breathing kagome lattices. Without loss of generality, we use magnetic vortices to demonstrate the principle and as a proof of the concept. The collective motion of vortices is described by the generalized Thiele's equation including an inertial term and a non-Newtonian gyroscopic term. We compute the vortex gyration spectra and find that the system is nontrivial and supports the topological corner states when d 2 /d 1 1.2 (d 2 /d 1 1.2 or d 1 /d 2 1.2) for triangle-shape (parallelogram-shape) lattice. The non-topological Tamm-Shockley edge state is also observed. Here d 1 and d 2 are the distances between two kinds of nearest-neighbor vortices, as shown in Fig. 1(a) and Fig. 4(a) for triangle and parallelogram structures, respectively. We show that the topological corner states emerge near the gyration frequency of a single vortex, and are robust against structure defects and disorders. Micromagnetic simulations confirm the theoretical results. From an experimental point of view, magnetic soliton lattices can be easily fabricated by electron-beam lithography [43], compared with the complex fabrication processes of phononic and photonic crystals. Our findings shall encourage experimentalists to observe the predicted higher-order topological solitonic states.

II. GENERALIZED THIELE'S EQUATION
We consider a kagome lattice of magnetic nanodisks with vortex states. Figure 1(a) plots the lattice structure with alternate distance parameters d 1 and d 2 . To accurately obtain the gyration spectrum of the vortex array, we start with the generalized Thiele's equation [48]: where U j = R j − R 0 j is the displacement of the vortex core from the equilibrium position R 0 j , G = −4πQwM s /γ is the gyroscopic constant with Q = 1 4π dxdym · ( ∂m ∂x × ∂m ∂y ) being the topological charge [Q = −1/2 for the vortex configuration shown in Fig. 1(b)], m is the unit vector of magnetization, w is the thickness of nanodisk, M s is the saturation magnetization, γ is the gyromagnetic ratio, M is the effective mass of the magnetic vortex [50][51][52], and G 3 is the third-order gyroscopic coefficient [53][54][55]. The conservative force can be expressed as F j = −∂W/∂U j where W is the potential energy as a function of the vortex displacement: [49,56,57]. Here K is the spring constant which is determined by the identity ω 0 = K/|G|, ω 0 is the gyrotropic frequency of a single vortex (see Supplementary Note 1), I and I ⊥ are the longitudinal and transverse coupling constants, respectively. Imposing U j = (u j , v j ) and defining ψ j = u j + iv j , we have:D where the differential operatorD = iω 3 , and ξ l = (I l + I l ⊥ )/2G, with l = 1 (or l = 2) representing the distance d 1 (or d 2 ) between nearest neighbor vortices. It is worth mentioning that parameters I l and I l ⊥ can be obtained by calculating the eigenfrequencies of a two-vortex system with different combinations of vortex polarities (see Supplementary Note 2). θ jk is the angle of the directionê jk from x-axis witĥ e jk = (R 0 k − R 0 j )/|R 0 k − R 0 j | and j is the set of all intracell and intercell nearest neighbors of j.

III. CORNER STATES AND PHASE DIAGRAM
It is known that the coupling strength I and I ⊥ strongly depends on the parameter d (d = d /r with d the distance between two vortices and r being the radius of nanodisk) [58][59][60]. Consequently, to calculate the spectrum and the phase diagram of vortex gyrations, we need to know the analytical expression of I (d) and I ⊥ (d). With the help of micromagnetic simulations for a simple system consisting of two nanodisks, we obtain the best fit of the numerical data:  Fig. 2(a) with symbols and curves representing simulation results and analytical formulas, respectively. In the calculations, we have adopted the material parameters of Permalloy (Py: Ni 80 Fe 20 ) [61,62] with G = 3.0725 × 10 −13 Js/m 2 . The spring constant K, mass M, and non-Newtonian gyration G 3 are obtained by analyzing the dynamics of a single vortex confined in the nanodisk [54,63]: K = 1.8128 × 10 −3 J/m 2 , M = 9.1224 × 10 −25 kg, and G 3 = 4.5571 × 10 −35 Js 3 /m 2 (see Supplementary Note 1). Then, by solving Eq. (2) numerically, we obtain the eigenfrequencies of vortex gyrations in the breathing kagome lattice. Figure 2(b) shows the eigenfrequencies of the triangle-shape system for different values d 2 /d 1 with a fixed d 1 = 2.2r. By studying the spatial distribution of the corresponding eigenfunctions, we find that corner states can exist only if d 2 /d 1 1.2, as indicated by the red line segment. Different choices of d 1 gives almost the same conclusion (see Supplementary Note 3). Furthermore, we calculate the phase diagram by systematically changing d 1 and d 2 . It is shown that the boundary separating topologically non-trivial and metallic phases lies in d 2 /d 1 = 1.2, while topologically trivial and metallic phases are separated by d 1 /d 2 = 1.2, as shown in Fig. 2(c). When d 2 /d 1 1.2, the system is topologically non-trivial and can support second-order topological corner states. The system is trivial without any topological edge modes if d 1 /d 2 1.2. Here, the trivial phase is the gapped (insulator) state, the metallic/conducting phase represents the gapless state such that vortices oscillations can propagate in the bulk lattice, and the non-trivial phase means the second-order corner state surviving in a gapped bulk. It is worth mentioning that the critical condition (d 2 /d 1 1.2) for HOTIs may vary with respect to materials parameters. For example, the critical value will slightly increase (decrease) if the radius of the nanodisk increases (decreases).
Topological corner states should be robust against disorders in the bulk but sensitive to them at corners. To verify these properties of corners states in our system, we calculate the eigenfrequencies of the triangle-shape breathing kagome lattice of vortices under bulk disorders of different strengths (Disorders at three corners are discussed in Supplementary Note 4), as shown in Fig. 2(d), where d 1 = 2.08r and d 2 = 3.60r (d 2 /d 1 = 1.73 > 1.2). Here, disorders are introduced by assuming the resonant frequency ω 0 with a random shift, i.e., ω 0 → ω 0 + δZω 0 , where δ indicates the strength of the disorder and Z is a uniformly distributed random number between −1 to 1. We average the calculation after 100 realizations of uniformly distributed disorders. Gaussian distribution of Z leads to similar results. We can see from Fig. 2(d) that with the increasing of the disorder strength, the spectrum of both edge and bulk states is significantly modified, while the corner states are quite robust. Furthermore, the artifact effect of the corner states are also discussed in Supplementary Note 4. These findings echo the observations in photonic and phononic systems [13][14][15][16][17][18][19][20][21].
We choose the same geometric parameters as . It is found that there are three degenerate modes with the frequency 927.6 MHz, represented by red balls. Then, we confirm that these modes are indeed secondorder topological states (corner states) by showing the spatial distribution of vortex gyrations in Fig. 3(d) with oscillations being highly localized at the three corners. Besides these findings, we also identify the edge states, denoted by blue balls in Fig. 3(a). The spatial distribution of edge oscillations are confined on three edges, as shown in Fig. 3(c). However, these edge modes are Tamm-Shockley type [64,65], not chiral. They propagate in both directions, that is confirmed in micromagnetic simulations (see Supplementary Note 5). Bulk modes are plotted in Figs. 3(b) and 3(e), where corners do not participate in the oscillations.
The higher-order topological properties can be interpreted in terms of the bulk topological index, i.e., the polarization [66,67]: where S is the area of the first Brillouin zone, A j = −i ψ|∂k j |ψ is Berry connection with j = x, y, and ψ is the wave function for the lowest band. We have numerically calculated the polarization and find (P x , P y ) = (0.499, 0.288) for d 1 = 2.08r and d 2 = 3.60r and (P x , P y ) = (0.032, 0.047) for d 1 = 3r and d 2 = 2.1r. The former corresponds to the topological insulating phase while the latter is for the trivial phase. Theoretically, the polarization (P x , P y ) is identical to the Wannier center, which is restricted to two positions for insulating phases. If Wannier center coincides with (0, 0), the system is in the trivial insulating phase and no topological edge states exist. Higher-order topological corner states emerge when the Wannier center lies at (1/2, 1/2 √ 3) [10,15]. The distribution of bulk topological index is consistent with the computed phase diagram Fig. 2(c).
For completeness, we also study the corner states in another type of breathing kagome lattice of vortices (parallelogramshape), with the sketch plotted in Fig. 4(a). We consider the same parameters as those in the triangle-shape lattice. It is interesting to note that the results of triangle-shape and parallelogram-shape lattices are closely related. Two opposite acuted-angle corners in the parallelogram are actually not equivalent: one via d 1 bonding while the other one via d 2 bonding; see Fig. 4(a). Only the d 2 bonding (bottomright) corner in the parallelogram-shape lattice is identical to three corners in the triangle-shape lattice. Therefore, for parallelogram-shape lattice, we can observe only one corner state either in the bottom-right corner [when

IV. MICROMAGNETIC SIMULATIONS
To verify the theoretical predictions of HOTIs above, we perform full micromagnetic simulations by considering a breathing kagome lattice consisting of a few identical magnetic nanodisks in vortex states, as shown in Fig. 1(a) and Fig. 4(a), with the same geometric parameters as those in Fig. 3 and Fig. 4, respectively. Micromagnetic software MUMAX3 [68] is used to simulate the dynamics of vortices in Py. The material parameters are as follows: the saturation magnetization M s = 0.86 × 10 6 A/m, the exchange stiffness A = 1.3 × 10 −11 J/m, and the Gilbert damping constant α = 10 −4 (in order to observe the vortex oscillations clearly, we have chosen a rather small damping parameter). In the simulations, we set the cell size to be 2 × 2 × 10 nm 3 . To excite the full spectrum (up to a cut-off frequency) of the vortex oscillations, we apply a sinc-function magnetic field along the x-direction with H 0 = 10 mT, f = 20 GHz, and t 0 = 1 ns, as plotted in Fig. 1(c). The exciting field is applied over the whole system. The spatiotemporal evolutions of the vortices center R j = (R j,x , R j,y ) in all nanodisks are recorded every 200 ps, with the total simulation time being 1000 ns. Here R j,x and R j,y are defined by R j,x = x|m z | 2 dxdy |m z | 2 dxdy and R j,y = y|m z | 2 dxdy |m z | 2 dxdy , with the integral region being confined in the j-th nanodisk.
To identify the energy band of higher-order topological edge states in triangle-shape lattice, we compute the temporal Fourier spectrum of the vortex oscillations at different positions [labeled with arabic numbers 1, 2 and 3, see Fig. 1(a)]. for the corner has a very strong peak, which does not happen for the edge and bulk. We therefore infer that this is the corner-state band with oscillations localized only at three corners. Similarly, one can identify the frequency range which allows the bulk and edge states, as shown by shaded area with different colors in Fig. 5(a). To visualize the spatial distribution of vortex oscillations for different modes, we choose four representative frequencies: 940 MHz for the corner state, 842 MHz for the edge state, and both 769 MHz and 959 MHz for bulk states, and then stimulate their dynamics by a sinusoidal field h(t) = h 0 sin(2π f t)x with h 0 = 0.1 mT applied to the whole system for 100 ns. We plot the spatial distribution of oscillation amplitude in Figs. 5(b)-5(e). One can clearly see the corner state in Fig. 5(d), which is in a good agreement with theoretical results shown in Fig. 3(d) (theoretically calculated corner state locates at 927.6 MHz). Spatial distribution of vortices motion for bulk and edge states are shown in Figs. 5(b) and 5(c), respectively. It is noted that vortices at three corners in Fig. 5(e) also oscillate with a sizable amplitude, which is somewhat quite unexpected for bulk states. We attribute this inconsistency to the strong coupling (or hybridization) between the bulk and corner modes, since their frequencies are very close to each other, as shown in Figs.  3(a) and 5(a).
Like the triangle-shaped case, we have identified the corner, edge, and bulk states by micromagnetic simulation in parallelogram-shaped lattice with the same sinusoidal exciting fields applied to the whole system. We compute the temporal Fourier spectrum of the vortex oscillations at different positions [denoted with arabic numbers 1, 2 and 3, see  Fig. 6(d) shows only one corner state at only one (bottom-right) acute angle, which is in a good agreement with theoretical results shown in Fig. 4(e). Spatial distribution of vortices motion for bulk and edge states are shown in Figs. 6(b) and 6(c), respectively. Interestingly, the hybridization between the bulk mode and corner mode occurs as well in parallelogram-shaped breathing kagome lattice, see Fig. 6(e).

V. CONCLUSIONS
We have investigated the higher-order topological insulator in triangle-shaped and parallelogram-shaped breathing kagome lattice of magnetic vortices. Phase diagram including various solitonic states was obtained theoretically. It was found that the second-order topological corner state emerge only under a critical geometric parameter. We interpreted these results by the bulk topological index. Micromagnetic simulations were performed to confirm all theoretical predictions. We envision the existence of higher-order topological solitonic insulators in other type of lattices (breathing honeycomb, for instance), which is an interesting issue for future study. Identifying higher-order topological magnon insulator is also an open question. We believe that the findings presented in this work shall encourage experimentalists to find higher-order topological states in magnetic systems, within current technology reach.

Supplementary Information
We provide here information about the determination of the spring constant K, the vortex mass M, and the non-Newtonian gyration coefficient G 3 in Supplementary Note 1, the calculation of coupling parameters I and I ⊥ between two vortices in Supplementary Note 2, and the eigenfrequencies of triangle-shape breathing kagome vortex lattice for different lattice constants in Supplementary Note 3. For triangle-shape lattice, the numerical demonstration of the robustness of corner states against disorders and defects are provided in Supplementary Note 4. The Tamm-Shockley edge state is shown in Supplementary Note 5. Finally the robustness of the corner states and the phase diagram for parallelogram-shape lattice are discussed in Supplementary Note 6.

Supplementary Note 1
We determine the spring constant K, vortex mass M, and non-Newtonian gyration G 3 through the following relations [54,63]: ω 0 = K/G, G 3ω 2 = G, G 3 (∆ω + ω 0 ) = M, 2ω = ω 1 + ω 2 , and ∆ω = |ω 2 − ω 1 |, where ω 0 is the frequency of the gyroscopic mode, ω 1 and ω 2 are the frequencies of the other two higher-order modes with opposite gyration handedness [54]. For Py, G = 3.0725 × 10 −13 Js/m 2 . We first numerically obtain the excitation spectrum of an isolated nanodisk with a vortex, with results plotted in Supplementary Fig. 7. We confirm that the lowest mode is an anti-clockwise gyrotropic mode and the other two high-frequency modes have opposite directions of gyration [69]. The frequency of the three peaks reads ω 0 = 2π × 0.939 GHz, ω 1 = 2π × 11.945 GHz, and ω 2 = 2π × 14.192 GHz. Substituting into the formula above, we get K = 1. To determine the two key parameters I and I ⊥ in the formula W = j KU 2 j /2 + j k U jk /2 with U jk = I U j U k − I ⊥ U ⊥ j U ⊥ k , U j = U j ·ê jk , and U ⊥ j = U j · (ẑ ×ê jk ), we consider a simplified system consisting of two nanodisks (with the same size as those in the main text), as shown in Supplementary Fig.  8. The eigenfrequencies of coupled two-vortex system can be expressed as ω = ω 0 (1 ± I /K)(1 ∓ P 1 P 2 I ⊥ /K) [58], where P 1 (or P 2 ) is either +1 or −1 depending on the vortex polarity. Therefore, once the frequencies of coupled modes for different combinations of vortex polarities (P 1 P 2 = 1 or P 1 P 2 = −1) are determined from micromagnetic simulation, we can derive I and I ⊥ according to the dispersion relation. Supplementary  Figs. 8

(a) and 8(b) [Supplementary Figs. 8(c) and 8(d)]
show the setup and spectra of the two-vortex system with d = 2.08r (d = 3.60r). By extracting the frequency corresponding to the peak and after some algebra, we obtain I 1 = 1.2894 × 10 −4 J/m 2 , I 1 ⊥ = 3.5849 × 10 −4 J/m 2 , I 2 = 2.1237 × 10 −5 J/m 2 , and I 2 ⊥ = 4.4399 × 10 −5 J/m 2 . Further, we have systematically calculated the d−dependence of coupling parameters I and I ⊥ , as shown in Fig. 2(a) in the main text. As mentioned in the main text, for triangle-shape breathing kagome lattice of vortices, the phases of the system are entirely determined by the ratio d 2 /d 1 . To further confirm this conclusion, we plot the eigenfrequencies of system for different value d 2 /d 1 , with four choices of d 1 , as shown in Supplementary Figs. 9(a)-9(d). It can be clearly seen that regardless of the value of d 1 , the system always shows corner states when d 2 /d 1 1.2. The effect of disorders in the bulk on eigenfrequencies of triangle-shape system have been disscussed in the main text. For comparison, Supplementary Fig. 10(c) shows the results when disorders are located at three corners, we can seen that the frequencies of the corner states have great changes. Besides, the larger the strength of disorders, the bigger the frequency difference of the three corner states. To have a better understanding on the role of disorder on corner states, we plot the eigenfrequencies of system for disorders in the bulk [Supplementary Fig. 10 fects. Then we apply a sinusoidal field h(t) = h 0 sin(2π f t)x with h 0 = 0.1 mT and f = 940 MHz to the whole system. Supplementary Figs. 11(b), 11(d), and 11(f) show the spatial distribution of oscillation amplitude with different defect types. Micromagnetic simulations are run for 100 ns. One can clearly see that the corner states are robust for those defects. Interestingly, some defects may introduce new corners and corner states, see Supplementary Fig. 11(d). But corner states always emerge at acuted-angle corners.

Supplementary Note 5
In previous sections, we have shown that there exist edge states along three boundaries. Here, we confirm that the observed edge modes are Tamm  mT, f = 842 MHz applied on one disk at the edge, indicated by the blue arrow in Supplementary Fig. 12, with the same parameters as those in the main text. Then we draw the spatial distribution of oscillation amplitude, it can be seen clearly that the propagation of vortex oscillations is bidirectional. This nonchiral mode can be simply explained in terms of the Tamm-Shockley mechanism [48,64,65].