Higher-order topological solitonic insulators

Abstract

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.

Introduction

Topological insulators are receiving considerable attention owing to their peculiar properties, such as chiral edge states, and potential applications in spintronics and quantum computing.^1,2,3,4 A conventional $n$-dimensional topological insulator only has $(n-1)$-dimensional (first-order) topological edge/surface modes according to the bulk-boundary correspondence.^1,2 However, a higher-order, e.g., kth-order, topological insulator allows $(n-k)$-dimensional topological boundary states with $2\le k\le n$, which goes beyond the standard bulk-boundary correspondence and is characterized by the bulk topological index.^{5,6,7,8,9,10,11,12} Interestingly, the experimental evidences of higher-order topological insulators (HOTIs) were reported so far only in classical mechanical and electromagnetic metamaterials.^{13,14,15,16,17,18,19,20,21} In terms of applications of HOTIs in spintronics, it is intriguing to ask if they can exist in magnetic system which is intrinsically, however nonlinear, in contrast to its phononic and photonic counterparts.

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}\ge 1.2$ (${d}_{2}/{d}_{1}\ge 1.2$ or ${d}_{1}/{d}_{2}\ge 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 Figs. 1a and 4a 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,⁴³ 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.

Fig. 1

Triangle-shape breathing kagome lattice of vortices. a Illustration of the triangle-shape breathing kagome lattice including 45 nanodisks of the vortex state. ${d}_{1}$ and ${d}_{2}$ are the distance between two nearest-neighbor vortices. Arabic numbers 1, 2 and 3 denote the positions of spectrum analysis for bulk, edge and corner states, respectively. b Zoomed-in details of a nanodisk with the radius $r=50$ nm and the thickness $w=10$ nm. c Time dependence of the sinc-function field $H(t)$ applied to the whole system

Results

Generalized Thiele’s equation

We consider a kagome lattice of magnetic nanodisks with vortex states. Figure 1a 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:⁴⁸

$${G}_{3}\hat{z}\times \frac{{{\mathrm {d}}}^{3}{{\bf{U}}}_{j}}{{\mathrm {d}}{t}^{3}}-M\frac{{{\mathrm {d}}}^{2}{{\bf{U}}}_{j}}{{\mathrm {d}}{t}^{2}}+G\hat{z}\times \frac{{\mathrm {d}}{{\bf{U}}}_{j}}{{\mathrm {d}}t}+{{\bf{F}}}_{j}=0,$$

(1)

where ${{\bf{U}}}_{j}={{\bf{R}}}_{j}-{{\bf{R}}}_{j}^{0}$ is the displacement of the vortex core from the equilibrium position ${{\bf{R}}}_{j}^{0},G=-4\pi Qw{M}_{{\mathrm {s}}}$/$\gamma$ is the gyroscopic constant with $Q=\frac{1}{4\pi } \iint {\mathrm {d}}x{\mathrm {d}}y{\bf{m}}\cdot (\frac{\partial {\bf{m}}}{\partial x}\times \frac{\partial {\bf{m}}}{\partial y})$ being the topological charge ($Q=+1/2$ for the vortex configuration shown in Fig. 1b), ${\bf{m}}$ is the unit vector of magnetization, $w$ is the thickness of nanodisk, ${M}_{{\mathrm {s}}}$ is the saturation magnetization, $\gamma$ 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 ${{\bf{F}}}_{j}=-\partial W/\partial {{\bf{U}}}_{j}$, where $W$ is the potential energy as a function of the vortex displacement: $W={\sum }_{j}K{{\bf{U}}}_{j}^{2}/2+{\sum }_{j\ne k}{U}_{jk}/2$ with ${U}_{jk}={I}_{\parallel }{U}_{j}^{\parallel }{U}_{k}^{\parallel }-{I}_{\perp }{U}_{j}^{\perp }{U}_{k}^{\perp }$.^49,56,57 Here $K$ is the spring constant which is determined by the identity ${\omega }_{0}=K/| G|$, ${\omega }_{0}$ is the gyrotropic frequency of a single vortex (see Supplementary Note 1), ${I}_{\parallel }$ and ${I}_{\perp }$ are the longitudinal and transverse coupling constants, respectively. Imposing ${{\bf{U}}}_{j}=({u}_{j},{v}_{j})$ and defining ${\psi }_{j}={u}_{j}+i{v}_{j}$, we have

$$\hat{{\mathcal{D}}}{\psi }_{j}={\omega }_{0}{\psi }_{j}+{\sum\limits_{k\in \langle j\rangle}}({\zeta }^{l}{\psi }_{k}+{\xi }^{l}{{\mathrm {e}}}^{i2{\theta }_{jk}}{\psi }_{k}^{* }),$$

(2)

where the differential operator $\hat{\mathcal{D}}=i{\omega }_{3}\frac{d^3}{{\mathrm {d}}{t}^{3}}-{\omega }_{M}\frac{{d^2}}{{\mathrm {d}}{t}^{2}}-i\frac{d}{{\mathrm {d}}t}$, ${\omega }_{3}={G}_{3}/\vert G \vert$, ${\omega }_{M}=M/\vert G \vert$, ${\zeta }^{l}=({I}_{\parallel }^{l}-{I}_{\perp }^{l})/2|G|$, and ${\xi }^{l}=({I}_{\parallel }^{l}+{I}_{\perp }^{l})/2|G$, 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}_{\parallel }^{l}$ and ${I}_{\perp }^{l}$ can be obtained by calculating the eigenfrequencies of a two-vortex system with different combinations of vortex polarities (see Supplementary Note 2). ${\theta }_{jk}$ is the angle of the direction ${\hat{e}}_{jk}$ from $x$-axis with ${\hat{e}}_{jk}=({{\bf{R}}}_{k}^{0}-{{\bf{R}}}_{j}^{0})/| {{\bf{R}}}_{k}^{0}-{{\bf{R}}}_{j}^{0}|$ and $\langle j\rangle$ is the set of all intracell and intercell nearest neighbors of $j$.

Corner states and phase diagram

It is known that the coupling strength ${I}_{\parallel }$ and ${I}_{\perp }$ strongly depends on the parameter $d$ ($d=d^{\prime} /r$ with $d^{\prime}$ 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}_{\parallel }(d)$ and ${I}_{\perp }(d)$. With the help of micromagnetic simulations for a simple system consisting of two nanodisks, we obtain the best fit of the numerical data: ${I}_{\parallel }=\mu_{0}{M}_{{\mathrm {s}}}^{2}r(-1.72064\times 1{0}^{-4}+4.13166\times 1{0}^{-2}/{d}^{3}-0.24639/{d}^{5}+1.21066/{d}^{7}-1.81836/{d}^{9})$ and ${I}_{\perp }=\mu_{0}{M}_{{\mathrm {s}}}^{2}r(5.43158\times 1{0}^{-4}-4.34685\times 1{0}^{-2}/{d}^{3}+1.23778/{d}^{5}-6.48907/{d}^{7}+13.6422/{d}^{9})$, as shown in Fig. 2a 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\times 1{0}^{-13}\ {{{{\rm{J}}\,{\rm{s}}\,\rm{rad}^{-1}\,{\rm{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\times 1{0}^{-3}$ J m${}^{-2}$, $M=9.1224\times 1{0}^{-25}$ kg, and ${G}_{3}=-4.5571\times 1{0}^{-35}$ J s³rad⁻³m⁻² (see Supplementary Note 1). Then, by solving Eq. (2) numerically, we obtain the eigenfrequencies of vortex gyrations in the breathing kagome lattice. Figure 2b 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}\ge 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. 2c. When ${d}_{2}/{d}_{1}\ge 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}\ge 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}\ge 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).

Fig. 2

Corner states and phase diagram for triangle-shape lattice. a Dependence of the coupling strength ${I}_{\parallel }$ and ${I}_{\perp }$ on the vortex–vortex distance $d$ (normalized by the disk radius $r$). Pentagrams and circles denote simulation results and solid curves represent the analytical fitting. b Eigenfrequencies of collective vortex gyration under different ratio ${d}_{2}/{d}_{1}$ with the red segment labeling the corner state phase. c The phase diagram. d Eigenfrequencies of the breathing kagome lattice of vortices under different disorder strengths

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. 2d, 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 ${\omega }_{0}$ with a random shift, i.e., ${\omega }_{0}\to {\omega }_{0}+\delta Z{\omega }_{0}$, where $\delta$ indicates the strength of the disorder and $Z$ is a uniformly distributed random number between $-1$ and $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. 2d 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 Fig. 2d to explicitly visualize the corner states and other modes in the phase diagram. In this case, from Fig. 2a, we have ${I}_{\parallel }^{1}=1.2894\times 1{0}^{-4}$ J m⁻², ${I}_{\perp }^{1}=3.5849\times 1{0}^{-4}$ J m⁻², ${I}_{\parallel }^{2}=2.1237\times 1{0}^{-5}$ J m⁻², and ${I}_{\perp }^{2}=4.4399\times 1{0}^{-5}$ J m⁻². The computed eigenfrequencies and eigenmodes are plotted in Fig. 3a, b–e. 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 second-order topological states (corner states) by showing the spatial distribution of vortex gyrations in Fig. 3d with oscillations being highly localized at the three corners. Besides these findings, we also identify the edge states, denoted by blue balls in Fig. 3a. The spatial distribution of edge oscillations are confined on three edges, as shown in Fig. 3c. 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 Fig. 3b, e, where corners do not participate in the oscillations.

Fig. 3

Eigenmodes in triangle-shape lattice. a Eigenfrequencies of triangle-shape kagome vortex lattice with ${d}_{1}=2.08r$ and ${d}_{2}=3.60r$. The spatial distribution of vortex gyrations for the bulk (b and e), edge (c), and corner (d) states

The higher-order topological properties can be interpreted in terms of the bulk topological index, i.e., the polarization:^66,67

$${P}_{j}=\frac{1}{S}\iint _{\text{BZ}}{A}_{j}{d}^{2}k,$$

(3)

where $S$ is the area of the first Brillouin zone, ${A}_{j}=-i\langle \psi | \partial {k}_{j}| \psi \rangle$ is Berry connection with $j=x,y$, and $\psi$ 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$\sqrt{3}$).^10,15 The distribution of bulk topological index is consistent with the computed phase diagram Fig. 2c.

For completeness, we also study the corner states in another type of breathing kagome lattice of vortices (parallelogram-shape), with the sketch plotted in Fig. 4a. We consider the same parameters as those in the triangle-shape lattice. Figure 4b shows the eigenfrequencies of system which are obtained by numerically solving Eq. (2). Interestingly, we see that there is only one corner state at the frequency equal to 927.6 MHz, represented by the red ball. Edge and bulk states are also observed, denoted by blue and black balls, respectively. To have a better understanding of these modes, we have plotted the spatial distribution of vortices oscillation, as shown in Fig. 4c–f. From Fig. 4e, one can clearly see that the oscillations for corner state are confined to one acute angle and the vortex at the position of two obtuse angles hardly oscillates. The spatial distribution of vortex gyration for edge and bulk states are plotted in Fig. 4c, d, and f. The robustness of the corner states and the phase diagram are discussed in Supplementary Note 6.

Fig. 4

Eigenmodes in parallelogram-shape lattice. a The sketch for parallelogram-shaped breathing kagome lattice of vortices. Arabic numbers 1, 2, and 3 denote the position of spectrum analysis for bulk, edge, and corner states, respectively. b Numerically computed eigenfrequencies for parallelogram-shaped system. The spatial distribution of vortices oscillation for the bulk (c and f), edge (d), and corner (e) states

It is interesting to note that the results of triangle-shape and parallelogram-shape lattices are closely related. Two opposite acute-angle corners in the parallelogram are actually not equivalent: one via ${d}_{1}$ bonding, while the other one via ${d}_{2}$ bonding (see Fig. 4a). Only the ${d}_{2}$ bonding (bottom-right) 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 ${d}_{2}/{d}_{1}\ge 1.2$; see Fig. 4f in the main text] or in the top-left corner [when ${d}_{1}/{d}_{2}\ge 1.2$; see Supplementary Fig. 7].

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 Figs. 1a and 4a, with the same geometric parameters as those in Figs. 3 and 4, respectively. Micromagnetic software MUMAX3⁶⁸ is used to simulate the dynamics of vortices in Py (see the “Methods” section).

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. 1a). Figure 5a shows the spectra, with black, blue, and red curves denoting the positions of bulk (Number 1), edge (Number 2), and corner (Number 3) bands, respectively. One can immediately see that, near the frequency of 940 MHz, the spectrum 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. 5a. 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 and 959 MHz for bulk states, and then stimulate their dynamics by a sinusoidal field ${\bf{h}}(t)={h}_{0}\sin (2\pi ft)\hat{x}$ with ${h}_{0}=0.1$ mT applied to the whole system for 100 ns. We plot the spatial distribution of oscillation amplitude in Fig. 5b–e. One can clearly see the corner state in Fig. 5d, which is in a good agreement with theoretical results shown in Fig. 3d (theoretically calculated corner state locates at 927.6 MHz). Spatial distribution of vortices motion for bulk and edge states are shown in Fig. 5b, c, respectively. It is noted that vortices at three corners in Fig. 5e 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. 3a and 5a.

Fig. 5

Micromagnetic simulation of excitations in triangle-shape structure. a The temporal Fourier spectrum of the vortex oscillations at different positions. The spatial distribution of oscillation amplitude under the exciting field of various frequencies, 769 MHz b, 842 MHz c, 940 MHz d, and 959 MHz e. Since the oscillation amplitudes of the vortex centers are too small, we have magnified them by 2 or 10 times labeled in each figure

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. 4a). The spectra are shown in Fig. 6a with the black, blue, and red curves indicating the positions of bulk (Number 1), edge (Number 2), and corner (Number 3) bands, respectively. Shaded area with different colors denote different modes. The spatial distribution of oscillation amplitude is plotted in Fig. 6b–e. Figure 6d shows only one corner state at only one (bottom-right) acute angle, which is in a good agreement with theoretical results shown in Fig. 4e. Spatial distribution of vortices motion for bulk and edge states are shown in Fig. 6b, c, respectively. Interestingly, the hybridization between the bulk mode and corner mode occurs, as well in parallelogram-shaped breathing kagome lattice (see Fig. 6e).

Fig. 6

Micromagnetic simulation of excitations in parallelogram-shape structure. a The temporal Fourier spectrum of the vortex oscillations at different positions. The spatial distribution of oscillation amplitude under the exciting field with different frequencies, 767 MHz b, 844 MHz c, 940 MHz d, and 964 MHz e. The simulation time is 100 ns. Since the oscillation amplitudes of the vortices centers are too small, we have magnified them by 2 or 10 times labeled in each figure

Discussion

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.

Methods

Numerical simulations

The numerical simulations are carried out by Mumax package. The material parameters are as follows: the saturation magnetization ${M}_{{\mathrm {s}}}=0.86\times 1{0}^{6}$ A m⁻¹, the exchange stiffness $A=1.3\times 1{0}^{-11}$ J m⁻¹, and the Gilbert damping constant $\alpha =1{0}^{-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\times 2\times 10$ nm³. To excite the full spectrum (up to a cut-off frequency) of the vortex oscillations, we apply a sinc-function magnetic field $H(t)={H}_{0}\sin [2\pi f(t-{t}_{0})]/[2\pi f(t-{t}_{0})]$ along the $x$-direction with ${H}_{0}=10$ mT, $f=20$ GHz, and ${t}_{0}=1$ ns, as plotted in Fig. 1c. The exciting field is applied over the whole system. The spatiotemporal evolutions of the vortices center ${{\bf{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}=\frac{\iint x| {m}_{z}{| }^{2}{\mathrm {d}}x{\mathrm {d}}y}{\iint | {m}_{z}{| }^{2}{\mathrm {d}}x{\mathrm {d}}y}$ and ${R}_{j,y}=\frac{\iint y| {m}_{z}{| }^{2}{\mathrm {d}}x{\mathrm {d}}y}{\iint | {m}_{z}{| }^{2}{\mathrm {d}}x{\mathrm {d}}y}$, with the integral region being confined in the $j$th nanodisk