Emergence of Meron Kekulé lattices in twisted Néel antiferromagnets

Abstract

A Kekulé lattice is an exotic, distorted lattice structure exhibiting alternating bond lengths, distinguished from naturally formed atomic crystals. Despite its evident applicability, the formation of a Kekulé lattice from topological solitons in magnetic systems has remained elusive. Here, we propose twisted bilayer easy-plane Néel antiferromagnets as a promising platform for achieving a “Meron Kekulé lattice”—a distorted topological soliton lattice comprised of antiferromagnetic merons as its lattice elements. We demonstrate that the cores of these merons are stabilized into the Kekulé-O pattern with different intracell and intercell bond lengths across moiré supercells, thereby forming a Meron Kekulé lattice. Moreover, the two bond lengths of the Meron Kekulé lattice can be fine-tuned by adjusting the twist angle and specifics of the interlayer exchange coupling, suggesting extensive control over the meron lattice configuration in contrast to conventional magnetic systems. These discoveries pave the way for exploring topological solitons with distinctive Kekulé attributes.

Introduction

A Kekulé lattice is a distorted lattice structure distinguished by alternating bond lengths. This lattice distortion has been suggested to give rise to intriguing physical phenomena not observed in undistorted lattices^{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27}. Notable examples include charge fractionalization^1,2, chiral-symmetry breaking^{3,4,5,6,7,8,9}, and emergence of flat bands¹² in graphene. Additionally, the appearance of topological band structures has been demonstrated in photonic^21,22 and phononic crystals^23,24. Experimentally, the Kekulé lattice has been realized in graphene^{9,11,25,26,27}. However, the concept of a Kekulé lattice can extend beyond these systems, holding promise for broader applications across various physical systems with distinct lattice elements. Particularly fascinating is the potential creation of a Kekulé lattice from topological solitons in magnetic systems, which would introduce a novel perspective, contrasting with the traditionally observed Bravais lattice forms of topological solitons, such as the triangular lattice seen in the Abrikosov vortex lattice in superconductors²⁸ and the skyrmion lattice in magnetic systems²⁹. This intriguing possibility has remained unexplored, representing a pivotal piece yet to be uncovered within the realms of Kekulé physics.

In this study, we discover the creation of a Kekulé lattice from antiferromagnetic merons in twisted bilayer easy-plane Néel antiferromagnets. Twist engineering of van der Waals materials represents an innovative approach to creating spatial modulation of interlayer coupling through moiré patterns, leading to novel emergent physical phenomena in two-dimensional materials^{30,31,32,33,34,35}. This approach enables the stabilization of topological solitons, including skyrmions and merons, in ferromagnetic systems^{35,36,37,38,39,40,41,42}. By conducting atomistic spin simulations on twisted Neél antiferromagnets, we demonstrate that the cores of the antiferromagnetic merons are stabilized due to the moiré-induced spatial modulation of interlayer exchange coupling⁴¹. Notably, the cores of these merons form a honeycomb-lattice-like structure across moiré supercells, exhibiting different intracell and intercell bond lengths (Fig. 1). This distorted structure deviates from commonly observed Bravais-lattice forms^{29,43,44,45,46}, yet aligning with the known characteristics of the Kekulé-O pattern⁹. Hence, we term our distorted meron lattice a “Meron Kekulé lattice." We illustrate that the two bond lengths of the Meron Kekulé lattice can be fine-tuned by manipulating the twist angle and specifics of interlayer exchange coupling, offering a high degree of control over the stable meron lattice configuration. These discoveries present a Meron Kekulé lattice as a new classification within Kekulé lattices, introducing a novel pathway to attain this structure. Moreover, the stability and adjustability of this distorted structure present intriguing possibilities for exploring topological solitons with precise control—a capability rarely achievable within conventional magnetic systems^{29,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59}.

Fig. 1: Schematic illustration of a Meron Kekulé lattice in a twisted bilayer Néel antiferromagnet.

Yellow dots denote the cores of antiferromagnetic (AFM) merons that form a Kekulé lattice structure. Black solid lines denote intracell bonds between meron cores within the same moiré supercell, while white solid lines denote intercell bonds between meron cores across different supercells. The white dashed line depicts a single moiré supercell. Blue color indicates parallel alignment, while red indicates antiparallel alignment between Néel vectors across the top and bottom layers. In the magnified image, arrows represent the in-plane components of Néel vectors, depicting the winding texture of a meron pair.

Results

AFM domain array

We construct twisted bilayer Néel antiferromagnets by rotating two magnetic layers in a honeycomb lattice with a relative twist angle (Fig. 2a). These twisted magnets can be described using a Heisenberg spin model given by

$$\begin{array}{rcl}H&=& J\sum\limits_{l=t,b}\sum\limits_{\langle i,j\rangle }{{\boldsymbol{S}}}_{i}^{l}\cdot {{\boldsymbol{S}}}_{j}^{l}+A\sum\limits_{l=t,b}\sum\limits_{i}{\left({{\boldsymbol{S}}}_{i}^{l}\cdot \hat{{\boldsymbol{z}}}\right)}^{2}\\ &&+\sum\limits_{i,j}{J}_{ij}^{\perp }{{\boldsymbol{S}}}_{i}^{t}\cdot {{\boldsymbol{S}}}_{j}^{b}.\end{array}$$

(1)

Here, ${{\boldsymbol{S}}}_{i}^{l}$ represents the spin at site i on the top layer (l = t) and the bottom layer (l = b). The parameter J = 1.0 meV represents the intralayer AFM exchange interactions between nearest-neighbor spins. A = 0.2 meV represents the single-ion anisotropy energy favoring in-plane spin orientations. ${J}_{ij}^{\perp }$ represents the interlayer AFM exchange interactions that modulate as a function of the coordinate displacements ${{\boldsymbol{r}}}_{ij}^{tb}={{\boldsymbol{r}}}_{i}^{t}-{{\boldsymbol{r}}}_{j}^{b}$ between two spins at i- and j- sites on the top and bottom layers, respectively. To describe the decaying behavior of ${J}_{ij}^{\perp }$ as a function of $| {{\boldsymbol{r}}}_{ij}^{tb}|$, we employ the following exponential function^{35,36,38,39,40}:

$${J}_{ij}^{\perp }={J}_{0}^{\perp }\exp \left[-\alpha \left(\left\vert {{\boldsymbol{r}}}_{ij}^{tb}\right\vert /d-1\right)\right].$$

(2)

Here, the parameter ${J}_{0}^{\perp }$ represents the maximum value of ${J}_{ij}^{\perp }$ at $| {{\boldsymbol{r}}}_{ij}^{tb}| =d$, where d denotes the perpendicular interlayer distance; the parameter α describes the decay rate of ${J}_{ij}^{\perp }$ with respect to the increase of $| {{\boldsymbol{r}}}_{ij}^{tb}|$. The twist angle θ defines the size of the moiré superlattice L through the relationship: $L\approx \frac{\sqrt{3}a}{\theta }$ (θ: in radians). We choose α = 15, ${J}_{0}^{\perp }=0.3\,{\rm{meV}}$, d = 7Ã…, and a = 4Ã…, aligned with typical values observed in diverse vdW magnetic materials⁶⁰. We use the twist angle θ = 1.61°. These parameter values are utilized throughout this study unless explicitly specified.

Fig. 2: Antiferromagnetic domain array.

a Moiré superlattice for a twist angle of θ = 5.08°. Colored circles highlight distinct local stacking patterns: AA (red), AB (blue), and BA (cyan). The yellow rhombus and black arrows denote the unit cell and lattice vectors of the moiré superlattice, respectively. b Schematic illustration depicting different local spin configurations (A-type and B-type Néel orders) preferred in each stacking order. Here, red and blue arrows represent spins on the A and B sublattice, respectively. Green arrows depict Neél vectors in each layer. c Interlayer exchange energy map (${J}_{i}^{\perp }$) within the B-type order for θ = 1.61°. Here, red and blue colors indicate the preference for the A-type and B-type orders, respectively. d Zero temperature magnetic phase diagram as a function of twist angle (θ) and intralayer nearest-neighbor exchange interaction (J), displaying AFM and magnetic domain (MD) phases. The order parameter Ψ_t is defined as ${\Psi }_{t}=\frac{1}{{N}_{t}}\left\vert {\sum }_{i}{{\boldsymbol{N}}}_{i}^{t}\right\vert$ (${{\boldsymbol{N}}}_{i}^{t}$: a Néel vector on each i-site, N_t: the number of sites on the top layer), where Ψ_t = 1 indicates the AFM phase, and Ψ_t < 1 signifies the MD phase. The markers depict the boundary between the two phases. The dashed line represents a fitting curve defined as J = 162/θ². e–g AFM domain array configuration of the MD phase for θ = 1.61° and J = 1meV. In e–f, the color scales denote the phase angles (Ï•_t,b) of the normalized Néel vectors ${{\boldsymbol{N}}}_{t,b}=(\cos {\phi }_{t,b},\sin {\phi }_{t,b},0)$ in the top and the bottom layers, respectively. The arrows illustrate the direction of the Néel vectors in the plane. In g, the relative orientation of the Néel vectors between the two layers (N_t â‹… N_b), where red (blue) represents the A-type (B-type) Néel order.

The moiré superlattice of the twisted bilayer encompasses various local stacking patterns, such as AA, AB, and BA, within its supercell (Fig. 2a). This stacking modulation leads to diverse spin alignments between the top and bottom layers³⁵. Two characteristic alignment patterns emerge:

$$\begin{array}{rcl}\,\text{A}-\text{type}\,{\rm{N}}{\acute{{\rm{e}}}}{\rm{el}}\,\text{order}\,:\,{{\boldsymbol{N}}}_{i}^{t}\cdot {{\boldsymbol{N}}}_{i}^{b}&=&-1,\\ \,\text{B}-\text{type}\,{\rm{N}}{\acute{{\rm{e}}}}{\rm{el}}\,\text{order}\,:\,{{\boldsymbol{N}}}_{i}^{t}\cdot {{\boldsymbol{N}}}_{i}^{b}&=&+1.\end{array}$$

(3)

Here, ${{\boldsymbol{N}}}_{i}^{t}$ and ${{\boldsymbol{N}}}_{i}^{b}$ represent the Néel vectors on the top and bottom layers, respectively, defined as:

$${{\boldsymbol{N}}}_{i}^{t,b}=\frac{{{\boldsymbol{S}}}_{i,A}^{t,b}-{{\boldsymbol{S}}}_{i,B}^{t,b}}{\left\vert {{\boldsymbol{S}}}_{i,A}^{t,b}-{{\boldsymbol{S}}}_{i,B}^{t,b}\right\vert }.$$

(4)

In the AA stacking regions, the A-type order predominates since spins on the same sublattice in the top and bottom layers (e.g., ${{\boldsymbol{S}}}_{i,A}^{t}$ and ${{\boldsymbol{S}}}_{i,A}^{b}$ for the sublattice A) are coupled through the AFM coupling ${J}_{ij}^{\perp }$ (the left panel of Fig. 2b). Conversely, the AB stacking regions favor the B-type order, as ${{\boldsymbol{S}}}_{i,A}^{t}$ and ${{\boldsymbol{S}}}_{i,B}^{b}$ are primarily coupled, while ${{\boldsymbol{S}}}_{i,B}^{t}$ and ${{\boldsymbol{S}}}_{i,A}^{b}$ remain effectively decoupled due to their large separations (the middle panel). The BA stacking regions prefer the B-type order with similar considerations (the right panel). The map of local interlayer exchange energy, computed as ${J}_{i}^{\perp }={\sum }_{j}{J}_{ij}^{\perp }{{\boldsymbol{S}}}_{i}^{t}\cdot {{\boldsymbol{S}}}_{j}^{b}$, where ${{\boldsymbol{S}}}_{i,A}^{t}={{\boldsymbol{S}}}_{i,A}^{b}=(1,0,0)$ and ${{\boldsymbol{S}}}_{i,B}^{t}={{\boldsymbol{S}}}_{i,B}^{b}=(-1,0,0)$ reveals three distinct local regions (Fig. 2c): (i) “AFM patches" preferring the A-type order, as indicated by their AFM coupling character between two Néel vectors (i.e., ${J}_{i}^{\perp } > 0$), (ii) “FM patches" preferring the B-type order, as indicated by their FM coupling character (i.e., ${J}_{i}^{\perp } < 0$), and (iii) a neutral intermediate region lacking any specific preferred order, exhibiting negligible coupling (i.e., ${J}_{i}^{\perp }\approx 0$).

Our atomistic spin simulations on Eq. (1), conducted using an iterative optimization method^39,40,41, reveal a zero temperature magnetic phase diagram shown in Fig. 2d. In the simulations, random initial configurations within a single moiré supercell were relaxed, applying periodic boundary conditions at the unit cell boundaries. For detailed methods, refer to Methods. This diagram displays two distinct magnetic phases: AFM and magnetic domain (MD) phases. The AFM phase displays a uniform spin configuration characterized by the B-type Néel order. In contrast, the MD phase exhibits a nonuniform configuration (Fig. 2e–g), where the A-type Néel order emerges within each AFM patch, forming circular-shaped domains, whereas the B-type Néel order persists outside these patches (Fig. 2g). We term this spin configuration an “AFM domain array"⁴¹. The AFM domain array configuration minimizes ${J}_{i}^{\perp }$ across both the AFM and FM patch regions. This energy reduction outweighs domain wall energy in a small twist angle regime θ < θ_c, where ${\theta }_{c}=C\sqrt{{\bar{J}}_{\perp }/J}$ and ${\bar{J}}_{\perp }=0.3$ meV and C = 23.24^39,40,41, resulting in the emergence of the MD phase as the ground state (Fig. 2d).

Meron Kekulé lattice

We demonstrate that the AFM domain array, depicted in Fig. 2e-g, can accommodate merons as topological defects around its domain boundary, leading to a Meron Kekulé lattice. Merons are vortexlike topological solitons carrying half-skyrmion numbers^44,57, typically found in easy-plane magnets^{44,47,48,49,50,51,52,53,54,55,56,57,58,59}. The local density profile for the skyrmion number is calculated as:

$$q(x,y)=\frac{1}{4\pi }{\boldsymbol{N}}(x,y)\cdot ({\partial }_{x}{\boldsymbol{N}}(x,y)\times {\partial }_{y}{\boldsymbol{N}}(x,y)),$$

(5)

where N(x, y) represents the normalized interpolation of the Néel vectors. Furthermore, in-plane swirling textures away from the cores of merons are characterized by vorticity. The local density profile for the vorticity is calculated as⁶¹:

$$v(x,y)=\frac{1}{\pi }\hat{{\boldsymbol{z}}}\cdot ({\partial }_{x}{\boldsymbol{N}}(x,y)\times {\partial }_{y}{\boldsymbol{N}}(x,y)).$$

(6)

Figure 3 presents a stable spin configuration obtained through the relaxation of a random initial configuration. On the top layer, the Néel vectors exhibit six-fold in-plane swirling textures within a supercell, which are labeled by M₁, M₂, M₃, ${\overline{{\rm{M}}}}_{1}$, ${\overline{{\rm{M}}}}_{2}$, and ${\overline{{\rm{M}}}}_{3}$ (Fig. 3a). The map for v(x, y), presented in Fig. 3c, displays positive and negative values for (M₁, M₂, M₃) and (${\overline{{\rm{M}}}}_{1}$, ${\overline{{\rm{M}}}}_{2}$, ${\overline{{\rm{M}}}}_{3}$), respectively, identifying the former as vortices and the latter as antivortices. Furthermore, the map for q(x, y), depicted in Fig. 3d, demonstrates nonvanishing skyrmion densities with sign alternation for each core. These densities are locally integrated around each core region: Q = ∫dxdyq(x, y), yielding skyrmion numbers Q ≈ −1/2 for the vortices and Q ≈ +1/2 for the antivortices (for detailed methods, refer to Methods). These half-skyrmion number characteristics identify (M₁, M₂, M₃) as merons and (${\overline{{\rm{M}}}}_{1}$, ${\overline{{\rm{M}}}}_{2}$, ${\overline{{\rm{M}}}}_{3}$) as antimerons^44,57.

Fig. 3: Meron Kekulé lattice (MK) configuration.

a Néel vectors on the top layer (N_t), showing three merons (M₁, M₂, M₃) and three antimerons (${\overline{{\rm{M}}}}_{1}$, ${\overline{{\rm{M}}}}_{2}$, ${\overline{{\rm{M}}}}_{3}$). b Néel vectors on the bottom layer (N_b), displaying almost homogeneous alignment. In a-b, the color scales denote the phase angles (Ï•_t,b) of the normalized Néel vectors ${{\boldsymbol{N}}}_{t,b}=(\sin {\theta }_{t,b}\cos {\phi }_{t,b},\sin {\theta }_{t,b}\sin {\phi }_{t,b},\cos {\theta }_{t,b})$ in the top and the bottom layers, respectively. The black color denotes the out-of-plane component in the direction perpendicular to the plane (i.e., N_t = (0, 0, 1)). The arrows depict the interpolation of the in-plane components of the Néel vectors. c Vortex density (v) map corresponding to a, with red (blue) color indicating a positive (negative) winding number. d Skyrmion density (q) map corresponding to a, with red (blue) color indicating a positive (negative) skyrmion number. e The relative orientation of the Néel vectors between the two layers (N_t â‹… N_b), with blue (red) indicating parallel (antiparallel) alignment. In each panel, the dashed line depicts a supercell.

The cores of the merons and antimerons form six vertices of a hexagon along the boundary of the AFM patch (Fig. 3a). We refer to this distinctive configuration as a “Meron hexad." Despite its intricate texture, the Meron hexad attains an AFM domain array configuration that minimizes interlayer exchange energy across the bulk region (Fig. 3e). Furthermore, the constituent meron cores of Meron hexads form a honeycomb-lattice-like pattern, where the distance between two nearest-neighbor meron cores from the same hexad is shorter than that between cores from adjacent Meron hexads (Fig. 3a). This distorted structure is referred to as a Kekulé lattice with a Kekulé-O pattern⁹. Thus, we term this distorted meron lattice a “Meron Kekulé lattice (MK)." The Kekulé structure of the MK state is characterized by two key length parameters: the intracell bond length within the same Meron hexad (l) and the intercell bond length across different Meron hexads (${l}^{{\prime} }=L-2l$). The distinction between l and ${l}^{{\prime} }$ indicates the deviation from the conventional honeycomb lattice, giving rise to an alternative Kekulé structure. To quantify this deviation, we introduce a dimensionless constant λ, referred to as the “Kekulé constant:"

$$\lambda =\frac{{l}^{{\prime} }-l}{L}=1-\frac{3l}{L},$$

(7)

where λ = 0 and λ > 0 indicate regular honeycomb and Kekulé lattices, respectively.

We illustrate how variations in the parameters of the spin model influence the values of l and ${l}^{{\prime} }$, as well as the AFM patch size (r). Figure 4b shows the evolution of l and ${l}^{{\prime} }$ as functions of L. We first observe a linear increase in the AFM patch size (r) with L, which is evident since the increase in L signifies expansions in all local stacking regions, including AFM patches. The increase in r correlates with an increase in l, as the minimization of interlayer exchange energy dictates l = r⁴¹. In a large L regime, this relationship holds precisely; however, it progressively deviates in a small L regime. This deviation from l = r contributes to the enhancement of Kekulé distortion, represented by λ (see the inset). The linear increase in ${l}^{{\prime} }$ directly follows the increases in L and l, as described by the relationship ${l}^{{\prime} }=L-2l$.

Fig. 4: Manipulation of the Kekulé structure through parameter adjustments.

a Schematic illustration of the Meron Kekulé lattice structure. Yellow dots denote the position of the meron cores. Black lines denote intracell bonds between meron cores within a Meron hexad, while light gray lines denote intercell bonds between meron cores across different Meron hexads. l and ${l}^{{\prime} }=l-2L$ denote the intracell and intercell bond lengths, respectively. r denotes the size of AFM patches marked by purple circles. b–d Manipulation of the Kekulé lattice structure through three different parameters: b the moiré superlattice size (L), c the decay rate parameter of interlayer exchange interactions ${J}_{ij}^{\perp }$ (α), and d the maximum strength of ${J}_{ij}^{\perp }$ (${J}_{0}^{\perp }$). In b, the inset shows the value of the Kekulé constant ($\lambda =\frac{{l}^{{\prime} }-l}{L}$). The dashed lines represent empirical fitting curves, described by r = 0.28L + 0.05, l = 0.31L − 1.04, and ${l}^{{\prime} }=0.38L+2.09$. In c, the insets display ${J}_{i}^{\perp }$ for α = 10 (left) and α = 30 (right), respectively. The same color scale used in Fig. 2c is applied in the insets. In d, the dashed line represents the fitting curve for the l values. e Effective energy function E(l) for the MK state, plotted in a normalized form: $F(l/L)=\frac{E(l)-\min [E(l)]}{\max [E(l)]-\min [E(l)]}$. Here, `min' and `max' indicate the minimum and maximum values, respectively. Different colors indicate distinct values of ${\log }_{10}(p)$, where p represents the dimensionless ratio between the Coulomb interactions and the moiré field, which represents the two effective interactions influencing merons in the effective model. The parameter r/L = 0.25 is used. f Illustration of the stabilization condition of the MK state as a function of L and ${J}_{0}^{\perp }$. The MK state is stable in the blue region while unstable in the devoid region. The markers represent simulation data points, while the dashed line is a fitting curve defined as ${J}_{0}^{\perp }=40/{L}^{2}$.

We find that an increase in α leads to the contraction of r, resulting in a decrease in l while adhering to the relationship l ≈ r. Meanwhile, ${l}^{{\prime} }$ increases according to ${l}^{{\prime} }=L-2l$ with a fixed L (Fig. 4c). The contraction of r is attributed to significant suppression of ${J}_{ij}^{\perp }$ for large values of $| {{\boldsymbol{r}}}_{ij}^{tb}|$ away from the center of an AFM patch (see the insets). Moreover, an increase in ${J}_{0}^{\perp }/J$ initially leads to an increase in l when ${J}_{0}^{\perp }/J$ is sufficiently small (Fig. 4d). However, l reaches its maximum value l_max ≈ r, and further changes in ${J}_{0}^{\perp }/J$ do not affect l. We also find that reducing the ratio d/a lowers l by decreasing r, similar to the effects of α. Furthermore, varying A/J does not alter l, while extreme values lead to a transition to other magnetic states. Further details regarding this analysis of d/a and A/J can be found in Supplementary Note 2.

Table 1 presents Meron Kekulé lattices at various length scales (l, ${l}^{{\prime} }$, and L). These structures can be created by manipulating the twist angle (θ)^32,62,63 or by modifying the parameters of interlayer exchange coupling (α, d/a, and ${J}_{0}^{\perp }/J$). The latter approach can be implemented using different antiferromagnetic materials or external controls such as gate voltage, intercalation, adsorption, strain, or surface modification⁶⁴. The resulting Kekulé lattices exhibit a wide range of distortion, with λ values ranging from 0.05 to 0.34, highlighting the adaptability of the Meron Kekulé lattice.

Table 1 Examples of Kekulé structures (characterized by four parameters $l,{l}^{{\prime} },L$, and λ) in the MK state with various parameters

We emphasize that the discrepancy in meron bond lengths, l and ${l}^{{\prime} }$, indicates an alternation in the attractive interactions between meron cores within each bond of the Meron Kekulé lattice. In a typical magnet, a pair of a meron and an antimeron experiences an attractive interaction due to their opposite vorticities. At sufficiently large separations, the interaction energy follows a logarithmic function of the distance between the pair⁶⁵, indicating that the strength of this interaction is inversely proportional to the distance. Consequently, the attractive interactions in the intracell and intercell bonds can be expressed as $\sim \frac{1}{l}$ and $\sim \frac{1}{{l}^{{\prime} }}$, respectively, which differ from each other since $l\ne {l}^{{\prime} }$ in the Meron Kekulé lattice. This bond-dependent interaction can be characterized by the ratio of the two interactions, defined as $\kappa \equiv \frac{{l}^{{\prime} }}{l}=\frac{1+2\lambda }{1-\lambda }$, which reaches κ ≈ 3.2 for the optimal distortion λ = 0.42 observed in our simulation, highlighting a significant difference between the interactions in the intracell and intercell bonds.

We attribute the intricate behavior of l observed in Fig. 4b–d to the interplay between two magnetic interactions influencing merons: (i) attractive Coulomb interactions between merons and antimerons with opposite vorticities^65,66, induced by the intralayer exchange interactions, and (ii) an effective moiré field generated by the bottom layer through interlayer exchange interactions⁴¹. To elucidate their interplay, we utilize an energy function derived from an effective continuum model:

$$\begin{array}{rcl}E(l)&=&{C}_{1}J\ln \left[\sqrt{\frac{2l}{3L}}\frac{{\left({L}^{6}-{l}^{6}\right)}^{2}\left({L}^{6}-64{l}^{6}\right)}{{L}^{6}{\left({L}^{6}+27{l}^{6}\right)}^{2}}\right]\\ &&+\frac{{C}_{2}{J}_{0}^{\perp }}{{a}^{2}}{\left(l-r\right)}^{2}.\end{array}$$

(8)

Here, the first term describes the Coulomb interactions ranging up to the nearest-neighbor supercells, and the second term represents the moiré field. The constants C₁ and C₂ are numerical factors on the order of unity. The relative strength of these two interactions is represented by the dimensionless ratio:

$$p=\frac{2{C}_{2}{J}_{0}^{\perp }{L}^{2}}{{C}_{1}J{a}^{2}}.$$

(9)

The value of l, where the MK state is stabilized, denoted as l_*, can be determined by minimizing E(l) with respect to l. When p is large, l_* ≈ r (e.g., l ≈ 0.25L for r = 0.25L in Fig. 4e), indicating that the moiré field term dominates. Conversely, a reduction in p enhances the role of the Coulomb interaction term, leading to a decrease in l_*, as shown by the shifts in the local minimum of E(l) in Fig. 4e. We assert that this competition plays a crucial role in the intricate behavior of l observed in Figs. 4b–d. In particular, our analysis demonstrates that the l_* value derived from E(l) captures the qualitative features of the observed l values in the simulation, supporting this assertion (the dashed lines in Fig. 4d). For the derivation of E(l) and the fitting method for l, refer to Supplementary Note 3.

The effective potential E(l) offers insights into the underlying mechanisms that support the stability of merons within the MK state. In conventional magnets, merons tend to be unstable due to attractive Coulomb interactions, which lead to pair annihilation during energy relaxation processes. However, in our twisted magnets, the moiré field mitigates this instability by imposing an energy penalty for the deviation of the meron cores from the AFM domain boundary, leading to the energetic stability of the MK state, shown by the convexity of E(l) (Fig. 4e). Furthermore, the effective model predicts that the stabilization of the MK state can be achieved by increasing either the interlayer exchange coupling (${J}_{0}^{\perp }$) or the superlattice size (L), which is confirmed in our numerical simulations (Fig. 4f).

Response to external magnetic fields

To investigate a potential means to distinguish the MD and MK states by bulk measurements, we analyze the magnetic response of the MD and MK states to external magnetic fields, incorporating the Zeeman term defined as:

$${H}_{{\rm{Zeeman}}}=-\sum\limits_{l=t,b}\sum\limits_{i}{\boldsymbol{B}}\cdot {{\boldsymbol{S}}}_{i}^{l},$$

(10)

where ${\boldsymbol{B}}=B\hat{{\boldsymbol{z}}}$ denotes external magnetic fields in the out-of-plane magnetization. The initial MD and MK states, presented in Figs. 3 and 4, are relaxed under external magnetic fields with continuously varying strength B. From the relaxed spin configurations for each state, the magnetic susceptibility, χ ≡ dM/dB, is calculated as a function of B, where $M=\frac{1}{N}{\sum }_{i,l}{{\boldsymbol{S}}}_{i}^{l}\cdot \hat{{\boldsymbol{z}}}$ is the average magnetization in the out-of-plane magnetization. Figure 5a, b presents χ of the MD state. Near B = 0, χ initially decreases as ∣B∣ increases, while it shows an upturn in ∣B∣ > 2.22 meV, exhibiting shoulder-like features around ∣B∣ = 4.88 meV (Fig. 5a). Beyond these shoulders, χ sharply drops, and subsequently vanishes at ∣B∣ > B_sat. = 6.86 meV (Fig. 5b), where M saturates to ∣M∣ = 1.

Fig. 5: Magnetic susceptibility of the MD and MK states.

a Magnetic susceptibility (χ = dM/dB) of the MD state with respect to the out-of-plane magnetization (M) as a function of an external magnetic field in the out-of-plane direction (B). b Zoom-out view of a. In a, b, the yellow color highlights the persistence of the MD state. c Magnetic susceptibility (χ) of the MK state. Different colors indicate three different magnetic states: MK state (red), vortex state (blue), and MD state (yellow). d Zoom-out view of c.

Figure 5c and d presents χ of the MK state. In the low-field regime (∣B∣ < B_c1 = 0.52 meV), χ exhibits an increasing tendency with ∣B∣ in contrast to the MD state. This behavior is attributed to the canted spin configurations in the core region that enhance the magnetization induced by the field. However, entering into the intermediate-field regime (B_c1≤∣B∣ < B_c2 = 5.62 meV), χ undergoes a sharp drop and subsequently exhibits a decreasing tendency. Our findings reveal that at the threshold ∣B∣ = B_c1, the spin flipping occurs in the perpendicular direction for spins with out-of-plane components opposite to the field direction. The resulting state maintains in-plane swirling textures while lacking out-of-plane components in their Néel vectors. Consequently, this state retains only vorticities but loses half-skyrmion numbers characteristics, thus being named a “vortex state" (for the detailed analysis of this state, refer to Supplementary Note 5). Finally, upon entering the high-field regime (∣B∣≥B_c2), the vortex state transitions into the MD state through the annihilation of meron-antimeron pairs. At this regime, χ aligns with that of the MD state. We highlight the two anomalous discontinuities at the transition points (∣B∣ = B_c1 and ∣B∣ = B_c2) as discernible signatures of the formation of merons in the MK state.

Metastability

Despite our specific focus on the configuration depicted in Fig. 3, the MK state can manifest in a diverse array of 2⁶ potential configurations, each defined by two possible vorticities (p = ±1) for the six meron cores within a moiré supercell. Among these configurations, only 2⁵ could exhibit variations in energy levels since global polarity reversal has no impact on the energy state. Our analysis indicates that despite the multitude of configurations, they maintain almost identical energy levels, deviating by less than 10^−6 percentage from each other. Additionally, completely interchanging the positions of vortices and antivortices introduces a factor of 2, while swapping the top and bottom layers introduces another factor of 2 (note that the original MK state, with merons in the top layer, spontaneously breaks the top-bottom layer symmetry), resulting in a total of G = 2⁸ potential variants for the MK state, all possessing the same energy level. This suggests that each supercell can assume one of the G configurations with nearly identical energies. Consequently, a system comprising N_s supercells can present a total of ${G}^{{N}_{s}}$ distinct configurations, underscoring the vast diversity and potential intricacy of the MK state. For a more detailed analysis, refer to Supplementary Note 4.

The configurational diversity of the MK state suggests that it might be more favorable at finite temperatures compared to the ground MD state. To evaluate this, we calculate the expected occupation number of Meron hexads in each moiré supercell, denoted by f_MH, at temperature T. Considering two possibilities, the MD and MK configurations for each supercell, we express f_MH as: ${f}_{{\rm{MH}}}\simeq \frac{G{e}^{-{E}_{{\rm{MK}}}/{k}_{B}T}}{{e}^{-{E}_{{\rm{MD}}}/{k}_{B}T}+G{e}^{-{E}_{{\rm{MK}}}/{k}_{B}T}}$, where E_MK and E_MD represent the energies of the MK and MD states from a single supercell, respectively. In this calculation, we assume the energy of each supercell is independent of the presence or absence of a Meron hexad (an approximation that might hold unless a significant number of Meron hexads occupy the supercells). By utilizing the calculated value, E_MK−E_MD = 21.44 meV at θ = 1.02°, we determine f_MH = 0.32 at T = 40 K. Similarly, for E_MK − E_MD = 9.29 meV at θ = 0.50°, estimated from the fitting, we derive f_MH = 0.49 at T = 20 K (for further details, refer to Supplementary Note 1). These results imply that achieving the MK state at a finite temperature is feasible by either increasing the temperature or adjusting the twist angle at a small value⁶⁷.

Discussion

We propose that our theory can be applied to CoPS₃^68,69, which exhibits two essential factors for achieving the Meron Kekulé lattice: AFM intralayer exchange interactions and easy-plane magnetic anisotropy. Another potential application lies in XPS₃ (X = Mn, Ni, Fe) and MnPSe₃, where their intrinsic easy-axis magnetic anisotropy can be potentially altered into easy-plane anisotropy through experimental control, such as strain, gate voltage, or surface adsorption^{70,71,72,73,74}.

For experimental observations, we propose utilizing scanning magnetometry techniques with nitrogen-vacancy centers, as well as Lorentz transmission electron microscopy and magnetic transmission soft X-ray microscopy, to directly observe meron pairs ranging in size from 40 to 10 nm at different twist angles (θ = 0. 5° to 2°). Indirect measurements can involve detecting the anomalous multiple superlattice periodicity manifested by different bond lengths of merons through neutron scattering experiments. Additionally, anomalous kinks in the magnetization curve, as observed in Fig. 5, can offer valuable insights for such indirect measurements, serving as an indication of the presence of merons.

An intriguing research question to explore is how the bond-dependent interactions of merons affect phonon-like collective excitation modes associated with their core vibrations. One potential outcome is the formation of gapped Dirac energy spectra, specifically driven by hopping amplitude modulations in an effective tight-binding model, akin to observations in graphene^9,11,25. These phenomena may even exhibit topological band structures, as outlined in studies of photonic and phononic crystals^21,22,23,24. Additionally, merons are influenced by distinctive physical effects due to their topological nature, such as gyrotropic forces, which can also lead to the emergence of topological band structures in core vibrational modes^75,76. Consequently, the interplay between lattice distortion and the topology of spin textures presents a compelling avenue to expand Kekulé physics beyond traditional systems composed of featureless lattice elements like atoms.

Methods

Relaxation of spin configurations

To achieve relaxed spin configurations, we employed an “iterative optimization method"^39,40,41,42. We initially generated random spin configurations represented by normalized vectors n_i = S_i/∣S_i∣ for a classical spin vector S_i at each i-th site. An effective magnetic field h_i, acting on each S_i, was computed as ${{\boldsymbol{h}}}_{i}\equiv \frac{\partial H[\{{{\boldsymbol{n}}}_{i}\}]}{\partial {{\boldsymbol{n}}}_{i}}$. Subsequently, we iteratively updated n_i by subtracting h_i from n_i, resulting in a spin configuration with reduced magnetic energy at each iteration. The update rule is defined as:

$${{\boldsymbol{n}}}_{i}^{(s+1)}=\frac{{{\boldsymbol{n}}}_{i}^{(s)}-c{{\boldsymbol{h}}}_{i}^{(s)}}{\left\vert {{\boldsymbol{n}}}_{i}^{(s)}-c{{\boldsymbol{h}}}_{i}^{(s)}\right\vert },$$

(11)

where (s) denotes the iteration step, and c controls the update rate. The iteration proceeds until ${{\boldsymbol{h}}}_{i}^{(s)}$ becomes parallel to ${{\boldsymbol{n}}}_{i}^{(s)}$, indicating that further updates will not alter the spin configuration and ensuring minimal magnetic energy. To guarantee numerical convergence, we employed the condition $\frac{1}{N}{\sum }_{i}| {{\boldsymbol{h}}}_{i}^{\perp }| < 1{0}^{-10}$, where ${{\boldsymbol{h}}}_{i}^{\perp }={{\boldsymbol{h}}}_{i}-{{\boldsymbol{n}}}_{i}({{\boldsymbol{n}}}_{i}\cdot {{\boldsymbol{h}}}_{i})$. This convergence criterion aligns with the stability condition in the Landau-Lifshitz-Gilbert equation governing the magnetic dynamics. Throughout the iterative updates, we applied periodic boundary conditions to maintain consistency across the boundaries of the moiré unit cell.

Determination of the phase diagram

To obtain the phase diagram presented in Fig. 2d, we used random initial configurations denoted as $\{{{\boldsymbol{n}}}_{i}^{(0)}\}=\{(\frac{{u}_{i}}{{n}_{i}},\frac{{v}_{i}}{{n}_{i}},\frac{{w}_{i}}{{n}_{i}})\}$, where u_i, v_i, and w_i are random numbers within the range of [−1, 1] and ${n}_{i}=\sqrt{{u}_{i}^{2}+{v}_{i}^{2}+{w}_{i}^{2}}$. These initial configurations were relaxed using an iterative optimization method to obtain the relaxed configuration {n_i}. The ground state was identified as the lowest energy configuration among {n_i} for each parameter set (θ, J). From the resulting ground states, the order parameters Ψ_l, defined as

$$\begin{array}{r}{\Psi }_{l}=\left\vert \frac{1}{{N}_{l}}\sum\limits_{i}{{\boldsymbol{N}}}_{i}^{l}\right\vert \end{array}$$

(12)

were calculated for each layer l = t, b. Here, N_l is the number of spins in each layer l. The AFM phase is characterized by Ψ_t = Ψ_b = 1, while the MD phase signifies Ψ_t < 1 and Ψ_b < 1. Ψ_t was presented in Fig. 2d. It was found that the difference between Ψ_t and Ψ_b is less than 10^−13.

Computation of the vorticity and skyrmion densities

To compute the vorticity and skyrmion density maps of the MK state, depicted in Fig. 3c, d, we performed interpolation of the Néel vectors representing the spin configuration from the original honeycomb lattice to a square lattice consisting of N_sq sites. Subsequently, the vorticity and skyrmion densities were calculated using these interpolated Néel vectors via Eqs. (6) and (5), respectively, with N_sq = 201². This method yielded consistent values for vorticity and skyrmion numbers, as demonstrated by obtaining V₁ = 0.9985 and Q₁ = −0.4926 for the meron M₁ when employing N_sq = 1, 001². These values closely align with the theoretical expectations of V₁ = 1 and Q₁ = − 1/2. The positions of the meron cores, as depicted in Fig. 4b–d, were identified based on the peak locations of these maps.

Determination of various Kekulé structures from different parameter sets

To determine l of the MK state from different parameter sets, as presented in Fig. 4b–d, we utilized various spin configurations characterized by different l values. These initial trial configurations were generated through the interpolation of the original MK state to attain diverse initial l/L values ranging from 0.2 to 1.2 with an interval of 0.02. Subsequently, for each parameter set, the initial configurations were relaxed through the iterative optimization method to obtain stable spin configurations that potentially possess different l. From the relaxed spin configurations, an energy profile of the MK state, E_simul.(l), was computed as a function of l. E_simul.(l) typically exhibits a convex behavior with a local minimum. The predominant MK states for each parameter set were selected based on the minimized energy configuration, and the l values were computed from these states. Similarly, the stabilization condition of the MK state, as presented in Fig. 4f, was determined from the energy profile; the presence (absence) of the local minimum in this profile indicates the stable (unstable) nature of the MK state in each phase space. Further details are available in the Supporting Information.

Predictions made from the effective energy function

To derive the fitting curve in Fig. 4d, we identified a local minimum l_* of E(l) in Eq. (8) by solving the equation $\frac{dE({l}_{* })}{dl}=0$ under the constraint $\frac{{d}^{2}E({l}_{* })}{d{l}^{2}} > 0$. In the calculations, the parameters J = 1, L = 24.60, and a = 0.4 were used, and r = 7.22 were adopted from the simulation. The fitting parameter C ≡ C₂/C₁ = 19.74 was employed. The fitting curve in Fig. 4f was derived from the formula ${J}_{0}^{\perp }=D/{L}^{2}$ where $D=\frac{pJ{a}^{2}C}{2}$. In the fitting, the parameters J = 1, a = 0.4, and p = 23.6 were used, where the value of p was determined by the threshold value of p for the stabilization of the MK state as determined from our numerical investigations (further details are available in the Supporting Information). The fitting parameter C was estimated as C = 21.33, which is consistent with the estimation of C = 19.74 from Fig. 4d.