Unidirectional spin-wave edge modes in magnonic crystal

One of the most outstanding features of 2D topological insulators is the presence of topological edge states within the bandgap of bulk materials.^1–4 Protected by the time-reversal symmetry, these edge states are robust against backscattering (i.e., they are unidirectional), which makes them a very promising candidate for the information carriers in future data processing devices. Topological edge states, insensitive to defects, were demonstrated theoretically and experimentally in many wave-hosting systems, e.g., electron,⁵ electromagnetic,⁶ acoustic,^7,8 or mechanical.⁹ 

Topological magnonics is an emerging topic in modern magnetism. Until now, it mainly exploits the analogies of topological waves considered in other wave hosting systems and is mainly limited to basic theoretic model systems. However, the field of magnonics provides a platform for studying new kinds of phenomena that are specific to spin waves (SWs), due to the intrinsic nonreciprocal and nonlinear properties of SWs^10,11 or the reconfigurability of magnetic states.¹² Therefore, the synthetic magnetic structures will play an important role in magnonics and will contribute to the development of topological waves in general.

Despite several theoretical demonstrations of topological magnons^13–26 or experimental proofs of a topological gap in 2D honeycomb ferromagnets,^27–29 the direct experimental demonstration of a magnonic unidirectional edge state is still missing. The main obstacles preventing experimental observations are the difficulty of nano-fabrication of complex geometry, low coupling between elements, difficulty controlling the magnetic state, low sensitivity of the in-plane dynamical field component, high density of modes when scaling up structures, or high damping in metallic ferromagnetic materials with Dzyaloshinskii–Moriya interaction.

A specific class of magnonic crystals (MCs) hosting the unidirectional edge states includes ferromagnetic MCs composed of dipolarly coupled nanoelements with a unit cell, magnetized in the closure domain state with fixed chirality.^14,15 A precursor model of this type of MC was presented in Ref. 13 as a square lattice with the four spins in the square unit cell. These four spins are located equidistantly around the center of the unit cell, with magnetizations circulating in-plane with the given chirality. Dipolar coupling between the four spins within a single unit cell gives rise to four azimuthal SW modes with s, p_x − ip_y, p_x + ip_y and $dx2−y2$-like symmetries. The eigen frequencies of both p-like azimuthal modes are degenerated. In the MC, the dipolar coupling of azimuthal modes from neighboring cells induces four SW bands.

By applying a perpendicular magnetic field, the spins get tilted out of the plane and the eigen frequencies of azimuthal modes are modified. The degeneracy of p-like azimuthal modes is lifted, and with increasing field, the frequency order of the four azimuthal modes can change. In the MC, this reordering of azimuthal mode frequencies can be accompanied by gap closing and successive band inversion signalizing the possible topological phase transition. Indeed, as shown in Ref. 13, with increasing magnetic field, the MC exhibits several phases of topological order manifested as nonzero Chern numbers in the SW bands. As expected, in the corresponding MC with finite width, the unidirectional edge states exist in the bandgap between the topologically nontrivial bands.

Nevertheless, the extension of the four-spin model to the nano/micro-scale system is challenging and is endowed with some of the above mentioned issues, e.g., difficulty to fabricate the sample, control the magnetic state, or high density of modes when scaling up structures. Even the dipolar coupling of modes remained elusive until recent experimental demonstrations.³⁰ 

Here, we present a numerical demonstration of the nanoscale realization of the topological properties shown for the four-spin model. The geometry of our MC potentially solves the experimental obstacles, and, therefore, our MC is a strong candidate for the experimental proof of the magnonic topological edge states. Moreover, we demonstrate that the properties of the edge states, such as direction of propagation, can be controlled with the strength of an external magnetic field.

The unit cell of the considered square lattice MC is composed of a thin ferromagnetic square with two lines cut along its diagonals [see Fig. 1(a)]. Therefore, the unit cell is composed of four elementary building blocks, right-angled isosceles triangles, see Fig. 1(b). The fabrication of this type of geometry requires removing only straight lines from the ferromagnetic film, and it can be realized with MFM scratching, local anodic oxidation,³¹ or a focused ion beam.^32–34 The closure domain magnetic state of the unit cell might nucleate spontaneously¹² or be induced with local MFM field,^35,36 all optical switching³⁷ or by combination of the in-plane field and variation in shape anisotropy of composing elements.³⁸ Furthermore, an efficient dipolar dynamical magnetic coupling between the elements is expected due to the large surface area of facing sides of elements and their small distance given by the thickness of the cut lines of the order of 10 nm.

The building blocks of the unit cell (we assumed Py triangles) considered in this manuscript had a thickness, t = 15 nm, and the base of the triangle was set to a = 150 nm, see Fig. 1(b). We assume an artificially low value of damping in spectral analysis (α = 0.0001) and a realistic value of Py’s damping in the demonstration of edge wave propagation (α = 0.005). Due to the shape anisotropy, the magnetic ground state of an isolated triangle is given by two quasiuniform magnetic states oriented along the base of the triangle (with +x and −x orientations). We can expect that in this geometry, the lowest-frequency spin wave modes will exhibit quantizations along the magnetization direction (along the x-axis). With an increase in the triangle size, higher-order quantizations along the x-axis are expected to exhibit a lower frequency than the quasi-uniform mode, due to the specificity of dipole interactions.³⁹ However, the modes with a node in the center of the triangle (odd modes along the x-axis) will have a relatively large area of excitation within the narrower part of the triangle. These modes will be more confined, as they are the modes that have the maximum amplitude in the center of the triangle (even modes along the x-axis). Therefore, the increase in frequencies due to confinement along the y-direction will be more pronounced for odd modes. Additionally, the magnetic state is deformed from the uniform ±x state, and spins tend to align with the sides of the triangle. In order to satisfy the continuous rotation of spins between two sides of the triangle at the triangle apex, a spin orientation along the x-axis is formed, generating a strong local demagnetizing field. This demagnetizing field decreases the magnitude of the effective field and, correspondingly, lowers the spin precession frequency at the triangle apex.

As a result, the lowest frequency quasi-uniform triangle excitation with strong localization at the triangle apex is found at 6.17 GHz for a = 150 nm, see Fig. 1(c) for the amplitude distribution. It is well below the frequency 10.15 GHz of the 1st quantized mode. Therefore, in the square formed of four Py triangles [Fig. 1(a)], we can consider the formation of four azimuthal modes from the lowest-frequency quasi-uniform excitations in triangles. These four azimuthal modes can be characterized by azimuthal quantum numbers, where mode with m = 0 is a uniform SW, mode with m = 2 resembles a standing wave, and modes with m = −1 and m = 1 correspond to the SWs traveling around the unit cell in clockwise (CW) and counterclockwise (CCW) directions, respectively. The frequencies of the three lowest azimuthal modes (with m = 2, 1, and −1) are well separated from the frequencies of azimuthal modes formed by higher-order excitations of triangle. Phase profiles for these three azimuthal modes are plotted in Figs. 1(e)–1(g). Properties of azimuthal modes in square unit cells exhibit qualitatively the same behavior as azimuthal modes in a simple four-spin model,¹³ assuming the coupling between the triangles within unit cells is stronger than the coupling between the triangles from neighboring unit cells. The width of the square diagonal line was set to S_in = 6 nm, and the separation between the square units was much larger (S_out = 30 nm) to satisfy this requirement.

The frequency dependence as a function of the triangle base, a, of the three lowest-frequency azimuthal modes in an isolated square unit cell is presented in Fig. 1(d). The lowest-frequency mode is m = 2 [Fig. 1(e)], separated by a frequency difference Δf₂₁ from two degenerated modes with m = 1 [Fig. 1(f)] and m = −1 [Fig. 1(g)]. The value of Δf₂₁ is about 1 GHz, which is a clear manifestation of strong dynamical dipolar coupling between the triangles in a unit cell. Moreover, as Fig. 1(d) shows, this Δf₂₁ decreases by only about 40% when a is tripled from 100 to 300 nm. This is an important result that demonstrates the possibility for scaling up the triangles’ size without introducing higher-order quantization modes of the individual triangular elements, which could lead to overlapping between magnonic bands and lowering the width of the full topological gap, as discussed in the following part. In this manuscript, we limit our investigation to a = 150 nm in order to demonstrate the effect and perform large-scale simulations of finite MC.

When creating a MC based on the unit cell described in Sec. II A, the azimuthal modes between neighboring cells are coupled, giving rise to the spin wave bands. The calculated bandstructure of infinite 2D MC is presented in Fig. 2(a), where three bands originating from three azimuthal modes with m = 2, 1, and −1 are shown. The bandwidths (a range of frequencies within a band) of these bands are about 0.7 GHz, which confirms significant dynamical dipolar coupling between the triangles of neighboring cells. The lowest band corresponding to m = 2 is separated by a gap of 0.6 GHz from the partially overlapping bands given by m = ±1, which is consistent with the fact that frequencies of azimuthal modes m = ±1 are degenerated. In order to show that the two lowest bands stem from the azimuthal modes m = 2 and 1, we have analyzed the spatial profiles of corresponding Bloch waves at k = Γ and M. Red points and green points in Figs. 2(a)–2(c) mark the Bloch waves with spatial profiles the same as azimuthal modes m = 2 and m = 1, respectively, depicted in Fig. 1(e) [Fig. 1(f)].

Motivated by the four-spin model, we study the evolution of the bandstructure of MC in the out-of-plane magnetic field. Particularly, we are searching for the bandgap closings, which are typical indicators of band inversions and successive topological phase transitions in periodic crystals.⁴⁰ Figure 2(b) shows an example of the bands in the vicinity of band closing occurring around B_z = 150 mT, where bands with m = 2 and m = 1 almost touch at the Γ point. At a slightly higher magnetic field, the gap closes and reopens again [see Fig. 2(c)], but the bands are now topologically nontrivial, i.e., they are no longer constructed solely by a single azimuthal mode. Indeed, red and green points reveal that in the vicinity of Γ point, the dominant azimuthal modes defining the bands are inversed when compared to the topologically trivial case in Figs. 2(a) and 2(b). Now, both lowest bands are composed of a mixture of azimuthal modes m = 2, 1 in a ratio depending on the position in a Brillouin zone.

To confirm the topological character of the formed bands, we perform the analysis of the infinite MC (with 32 × 32 unit cells and periodic boundary conditions), extract the Bloch wavefunctions and perform numerical calculation of the Chern numbers (see Sec. IV). As expected, we found that the Chern numbers for the two lowest bands at B_z = 0 T are zero, i.e., these bands are topologically trivial. However, the Chern numbers for the two lowest bands at B_z = 330 mT are C = 1 and −1 for the lower and upper bands, respectively, which confirms their topological nature.

In order to demonstrate the topological properties of the magnonic gap observed in Fig. 2(c), we perform calculations of the dispersion relation of the magnonic ribbon, consisting of an infinite number of square elements along the longitudinal (x) direction and eight square elements along the transversal (y) direction. The computational cell is set to 32 × 8 with periodic boundary conditions in the x-direction, so that 32 points of the dispersion can be plotted. The results of these calculations are presented in Fig. 3(a), where the dispersion along the x axis is calculated at the bottom edge of the finite ribbon (left image) and along the center (right image). Two bands, separated by the bandgap of the width ∼200 MHz are found at the center of the structure and the unidirectional wave mode connecting the bands was found at the edge of the structure. Furthermore, we perform simulation of the same infinite ribbon at k = 0 point (number of cells: 1 × 8) with high density mesh along the z-axis (size of the computational mesh is c_x, c_y, and c_z = 3 nm) to verify the influence of the non-uniform mesh along the z-axis. The intensity of the wave modes is plotted in Fig. 3(b) and shows the presence of strong edge localization at comparable frequencies.

Finally, we perform simulation of the edge state propagation in the finite structure without periodic boundary conditions (10 × 8 elements) with one defect, see Fig. 3(c). The sine signal was applied to one unit cell placed in the top right corner at frequency f = 5.65 GHz, as indicated by the white dotted lines on Fig. 3(a). The intensity of the dynamical excitation is plotted after 10 ns of the maximum of the sine signal and clearly demonstrates clockwise unidirectional propagation without scattering at the defect. The realistic value of Py damping was taken into account in the simulation of edge-state propagation (α = 0.005).⁴¹ 

According to the four-spin model, the gaps closing in the bandstructure of MC appear as a consequence of the frequency reordering of azimuthal modes. Namely, with an increasing magnetic field B_z, the degeneracy of the azimuthal modes m = ±1 is lifted. The reason is the following: The unit cell in the closure domain state has a zero component of magnetization m_z at B_z = 0. As the magnetic field increases, the m_z component of magnetization is induced. The influence of m_z on the splitting of the degeneracy of azimuthal modes was explained for exchange spin waves in isolated rings in terms of the geometric Berry phase acquired by an azimuthal spin wave.^42,43 The same argument can be applied for the system considered here, i.e., the square composed of triangles with azimuthal modes formed by dipole interactions between triangles. As a result, the frequency of azimuthal mode m = 1 decreases with increasing B_z in contrast to the frequency of m = −1, which increases. Moreover, with increasing B_z, the frequencies of the two lowest azimuthal modes, with m = 1 and m = 2, approach each other and eventually overlap, which leads to hybridization and gap closing between corresponding bands in MC. In the four-spin model, when increasing B_z, the first such gap closing appears at the Γ point and then the closings at the X and M points follow. Each of these gap closings is associated with a change in the topological phase of MC.

In the case of our MC, the azimuthal modes and bands show qualitatively similar behavior as described by the four-spin model. To extend the analysis of the bandstructure evolution in magnetic field, we show in Figs. 4(a)–4(c) the frequencies of magnonic bands at specific points k = Γ, K, and M respectively. As the magnetic field increases in the range from 0 to 300 mT, the frequencies of the bands formed by m = 1 and m = −1 azimuthal modes decrease and increase, respectively, i.e., the gap between those bands expands. The magnetic field in the considered range has a relatively low influence on the band formed by the m = 2 azimuthal mode. As a result, at B_z = 200 mT, the two lowest bands touch at the Γ point [see Fig. 4(a)], and the magnonic crystal undergoes a topological phase transition into a topologically non-trivial state. With further increase in magnetic field, the magnonic gap between the two lowest bands reopens, and the bands become inverted in the vicinity of Γ point, as discussed in section B.

Elevating the magnetic field above 300 mT, the frequencies of all bands steeply decrease. Figure 4(b) shows that at B_z = 430 mT, there is another gap closing between the two lowest bands at the X point, signaling a topological phase transition. Here, the Chern numbers of the two lowest bands are reversed, i.e., from C = 1 (lowest band) and −1 (2nd lowest band) in the range 200 < B_z < 430 mT to C = −1 and 1. This reversal of Chern numbers has direct consequence on the direction of topological edge states in finite MCs, as discussed in detail in Sec. II E. When increasing the magnetic field even more, the gap of about 150 MHz between the two lowest bands reopens again, which is visible in the inset of Fig. 4(b).

Further increase of the field to B_z = 760 mT shows the crossing of the m = 2 and m = 1 bands at the k = M point, and it is followed by a change of magnetic order at B_z = 782 mT (from a chiral state with a nonzero out of plane component to a state magnetized fully out of plane). Within the region between B_z = 760 mT and B_z = 782 mT, we observe complex band structures due to multi-mode overlapping [see inset in Fig. 4(c)].

Qualitatively, the same results and order of topological phases were predicted in the four-spin model in Ref. 13. The tight binding model presented there explains the main features of the dipole interaction that lead to the presence of the topological phases and inversion of the Chern numbers at the k = X point for the considered MC.

Furthermore, we present the band dispersions of two ribbons with a straight edge and a zigzag edge [Figs. 5(a) and 5(e) respectively], calculated at the bottom edge of the ribbons, at different magnetic fields in the range from 330 to 530 mT. These magnetic field ranges cover the transition between two topological phases occurring at 430 mT. Figures 5(b) and 5(f) show that the clockwise propagating edge state with positive group velocity ∝df/dk is present in the gap at B_z = 330 mT. With an increase in the magnetic field toward the transition at B_z = 430 mT, the magnonic gap width is decreasing [see Figs. 5(c) and 5(g)]. After the transition, at higher fields B_z = 530 mT [Figs. 5(d) and 5(h)], the edge state is observed in the gap with the opposite propagation direction (negative group velocity) as compared to the field at B_z = 330 mT. The change of edge state propagation direction is triggered by the change of the topological phase, accompanied by the reversal of the Chern numbers after the crossing of bands at the k = X point [Fig. 4(b)].

Finally, based on the change of topological phase induced by a magnetic field, we demonstrate the effect of a magnonic switch: a magnonic element that transfers the information in the direction determined by the value of the magnetic field and frequency of excitation. We performed the simulations of finite structure with the zigzag edge and plot the signal transmitted at time frame 10 ns after the maximum of the sin impulse for two magnetic fields B_z = 330 mT [see Fig. 6(a)] and B_z = 530 mT [see Fig. 6(b)]. The excitation frequencies were 5.91 and 4.96 GHz, respectively. As expected, the direction of the edge state is reversed at B_z = 530 mT as compared to the direction of propagation at the field B_z = 330 mT. The control of topological state direction creates a base for the development of reprogramable, immune to back-scattering magnonic applications, e.g., highly efficient magnonic diodes.

We have numerically demonstrated a topological magnonic crystal (MC) hosting unidirectional edge states that has potential for experimental realization. The MC proposed in this manuscript is a nanoscopic extension of the four-spin model,¹³ and it simultaneously addresses most of the issues (e.g., complex geometry, weak coupling between the elements, high density of modes leading to a small magnonic gap, …) preventing experimental demonstration of the topological edge states.

Namely, the geometry of our MC is very simple—it can be created by removing narrow straight lines from the thin Permalloy film. The unit cell of the MC consists of four isosceles right-angle triangles magnetized along their bases. Static magnetic moments of triangles circulate around the center of the unit cell with fixed chirality. The lowest-frequency spin wave eigenexcitation of a single triangle is well separated from the higher modes, and its amplitude is localized at the triangle apex. Therefore, the excitations of four triangles within a unit cell are strongly dipolarly coupled, creating four azimuthal modes with significant frequency differences. Dipolar coupling between the azimuthal modes from neighboring cells gives rise to four magnonic bands in the corresponding MC, separated by a wide gap of about several hundreds of MHz.

The optimization of the structure dimensions for an experimental realization of topological MC with nanoetching techniques is needed. This task should take into account the reduction of intercell dipolar coupling with scaling up of the structure dimensions, e.g., an increase in the lattice constant or separation between elements. However, this effect can be compensated by an increase in the structure thickness.

We have shown that by increasing the perpendicular magnetic field, the MC undergoes a couple of topological phase transitions, and the unidirectional edge states can be excited within the bulk bandgap in the MC with a finite width. Moreover, by tuning the magnetic field and the excitation frequency, the propagation direction of the edge state can be reversed, which could be utilized in a magnonic switch. We consider that this property will allow for designs of reprogrammable magnonic elements based on robust topological waves.

Throughout the present paper, in the numerical simulations, we use the following values of material parameters: saturation magnetization M_s = 0.8× 10⁶ A/m, exchange stiffness constant A = 1.3× 10⁻¹¹ J/m. The Gilbert damping parameter used in the dynamical simulations was α = 0.0001 for the calculation of mode frequencies and dispersion relations and α = 0.005 for the confirmation of unidirectional edge propagation within the gap. This set of parameter values corresponds to Py films. The mesh size was set to c_x = c_y = 3 nm and c_z = 15 nm for most of the simulations, with the exception of the high mesh density simulation, where the mesh was set to c_x = c_y = c_z = 3 nm.

We investigate the spectrum of spin excitations by performing finite difference time domain (FDTD) micromagnetic simulations using mumax³.⁴⁴ The calculations are based on the Landau–Lifshitz (LL) equation describing the dynamics of reduced magnetization m given by

where τ is the torque, which can be expressed as

Here, γ is the gyromagnetic ratio, α is a dimensionless damping parameter, and B_eff is the effective magnetic field, which is the sum of the external magnetic field B_ext, the magnetostatic field B_m, the isotropic Heisenberg exchange field B_ex (proportional to the exchange stiffness constant A),

Our simulations were performed as follows: An approximate closure domain state was assumed as the initial magnetization configuration in the unit cell, which was then relaxed to the minimum energy state. The stable magnetization configurations were excited by a low-amplitude variable magnetic field (b_x, b_y, b_z), the time dependence of which is represented by the continuous function

with the cut-off frequency f_max = 12.5 GHz and t₀ = 1 ns. The modulation amplitude of the sinc signal was I = 0.5 mT. For efficient excitation and detection of the azimuthal modes, a random spatial modulation, RN(x, y), of the external magnetic field was introduced in each cell. The RN is characterized by a normal distribution centered at 0 with a standard deviation equal to 1.

In order to confirm the unidirectional propagation of the edge mode, the sin signal was applied locally in the longitudinal direction, i.e., it has the form of [b_x(t), 0, 0] with

where I_sin = 0.05 mT, t₀ = 10 ns, and f is the excited frequency.

In order to evaluate the spectral properties of periodic structures, we follow the procedure described in Ref. 45. We applied a multidimensional Fourier transform to the dynamical signal of magnetization components δm_x,y,z(x, y, t), i.e., FT_β(f, k_x, k_y) = FFTN(δm_β(x, y, t)) with β = x, y, z. By summing the squares of magnitudes of functions FT_β, we receive the power spectral density function, P(f, k_x, k_y), that is periodic in k-space with periodicity defined by the reciprocal lattice vector $G=(2πΛ,2πΛ)$⁠, where Λ = a + S_out is a period of MC. For each point (k_x, k_y) in the Brillouin zone, we identify the peaks in P(f, k_x, k_y), which yield the frequencies f_n(k_x, k_y) of individual spin wave bands numbered by n.

In order to evaluate the spectral properties of ribbons (i.e., structures periodic along the x direction and finite along the y direction), we performed the Fourier transform of δm_x,y,z(x, y, t) in the x direction and time and obtained FT_β(y, f, k_x). The summation of the squared magnitudes of FT_β over β yielded the ribbon power spectral density P(y, f, k_x). The maxima of P define the ribbon local bandstructure f_n(y, k_x) at different transversal positions y, as shown in Figs. 3 and 5.

To explore the topological properties of spin wave bands, the usual procedure is to translate the undamped (α = 0) spin wave dynamics driven by LL Eq. (1) into a formalism analogous to the quantum systems driven by the Schrödinger equation.^13,14,46 Then, the Bloch formalism and standard procedures for the calculation of topological invariants can be applied. At first, the LL equation is linearized, considering that the dynamical part of the magnetization, δm(r), is small in comparison to the static part, m₀(r). Then, the magnetization at each position r_j is written on a basis where m₀(r_j) is oriented along the z-axis and δm(r_j) lies in an xy-plane, i.e., δm(r_j) = (δs_x(r_j), δs_y(r_j), 0).

Assuming the spin excitations in a periodic magnonic crystal, the Bloch theorem can be applied, i.e., $δsx(rj)=Xk(rj)eik.rj$ and $δsy(rj)=Yk(rj)eik.rj$⁠, where X_k(r_j) and Y_k(r_j) are time-dependent and periodic with the periodicity of the magnonic crystal.

Next, the magnonic creation m₊(r_j) = X_k(r_j) + iY_k(r_j) and annihilation m₋(r_j) = X_k(r_j) − iY_k(r_j) fields are introduced, which enables to rewrite the linearized LL equation for magnonic crystal in a Schrödinger-like form,

where H is Hermitian and η is a matrix with only diagonal elements being 1, −1, 1, −1, ….

Vector $ψ(k)=(m+(r1),m−(r1),m+(r2),m+(r2),…)T$ has 2N elements, with N being the number of coupled spins in a unit cell of magnonic crystal. When assuming ψ(k) = φ(k)e^−iωt, the stationary solutions of (6) obey the equation

In general, the Chern number of the n-th band is defined as

where the integral runs over the whole Brillouin zone and F₁₂(k) is a Berry curvature defined as

with A_μ(k) being the Berry connection,

In the case of spin wave eigenstates φ(k) obeying Eq. (7), the scalar product in (10) reads

The expression for Berry connection (10) is, however, gauge-dependent, and, in practice, this formula can be used only in the case when the analytical form of eigenstates φ with a fixed gauge is known. In our case, we obtain the eigenstates φ(k) numerically from the micromagnetic simulation (see the paragraph below), and the complex phase of each eigenstate is arbitrary. In addition, the eigenstates φ(k) are evaluated at discrete points across the Brillouin zone; therefore, the integral in (8) is replaced with the sum,

where the summands can be expressed in a gauge invariant form^47,48

with

In order to obtain the Bloch wave-functions of the spin wave bands from our micromagnetic calculations, we extract the values of Fourier transforms $FTβ(fn(kx,ky),kx+r2πΛ,ky+s2πΛ)$⁠, where r, s are integers and set 0 for other values of FT_β(f, k_x, k_y) functions. Finally, we performed the inverse Fourier transformation to obtain the Bloch wave function $δm′n,kx,ky(x,y)$ of the n-th band at k = (k_x, k_y). Functions $δm′n,kx,ky(x,y)$ can be readily transformed into the form of φ(k), and the Chern number (12) of the n-th band can be calculated.

The research has received funding from the Slovak Grant Agency APVV (Grant No. APVV-19-0311 (RSWFA)) and the VEGA agency (Grant No. VEGA-2/0177/21). This study was performed during the implementation of the project Building-up Center for advanced materials application of the Slovak Academy of Sciences, ITMS project code 313021T081 supported by Research and Innovation Operational Program funded by the ERDF (50%). The research has received funding from the National Science Center of Poland, Grant No. UMO-2018/30/Q/ST3/00416. The simulations were partially performed at the Poznan Supercomputing and Networking Center (Grant No. 398).

The authors have no conflicts to disclose.

J. Feilhauer: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). M. Zelent: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Zhiwang Zhang: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). J. Christensen: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). M. Mruczkiewicz: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.