Design of 2D skyrmionic metamaterials through controlled assembly

Introduction

Skyrmionic textures represent swirling topological spin configurations, typically stabilized by the Dzyaloshinskii-Moriya interaction (DMI). There are three primary skyrmionic textures: skyrmions^1,2,3,4,5,6, antiskyrmions^{7,8,9,10,11,12}, and skyrmionium^13,14,15, distinguished by their topological charges, morphology, and properties (such as the presence or not of the skyrmion Hall effect¹⁶). Moreover, recent studies have identified complex, higher-order skyrmionic textures carrying an absolute topological charge greater than 1^17,18,19. The emergent physical properties of skyrmionic textures show great potential in the design of next-generation sustainable spintronic applications, such as neuromorphic computing²⁰, race-track memory²¹, data encoding, and storage²².

Advances in understanding interactions among skyrmionic textures, such as attraction, repulsion, and annihilation, resonate with the well-established principles of ’Click chemistry’²³. ’Click chemistry’ is a process characterized by the synthesis of substances through the union of smaller units like molecules, with a focus on controlled assembly^{24,25,26,27,28,29,30,31,32,33,34}. Expanding this concept to ’click magnetism’, we introduce ’skyrmionic controlled assembly’, utilizing skyrmions, antiskyrmions, and skyrmioniums as elemental skyrmionic building blocks. The simulated controlled assembly protocol involves manually manipulating the position of these skyrmionic building blocks, or known textures, to induce various assembly processes that give origin to new textures with arbitrary topological charges. This methodology not only facilitates the integration of textures, but also enables the creation of metastable extended magnetic textures from these building blocks.

These extended textures are similar to a class of artificial materials, namely metamaterials^35,36, which are defined by their custom-engineered configurations and bespoke building blocks, e.g., nanoscale wrinkles and split-ring resonators. Metamaterials are extensively applied in nanostructure designs in the fields of optics and plasmonics to develop sustainable applications which offer capabilities beyond those of conventional materials, such as a negative refractive index³⁷. In this work, similar to the definition of metamaterials, we propose a skyrmionic metamaterial concept designed to accelerate the discovery of novel promising magnetic textures with intricate topologies for spintronics applications. We elaborate on this concept by detailing the creation of (i) lattice-like skyrmionic metamaterials, defined as spin-texture lattices that possess a topological charge greater than 1 within their primitive unit cells, as well as (ii) flake-like skyrmionic metamaterials, which are large-scale high-order antiskyrmions with an underlying lattice structure, and (iii) cell-like skyrmionic metamaterials that incorporate a nesting of various texture types. Each of these magnetic textures showcases a non-trivial topology and is constructed from fundamental skyrmionic units. Based on the aforementioned ’Click chemistry’ notion applied to skyrmionic textures, the key concept of our systematic exploration of skyrmionic metamaterials is also the controlled assembly, which enables the artificial placement of these elemental skyrmionic building blocks.

Employing spin dynamics simulation tools in our controlled assembly protocols, we investigate the relative stability of the skyrmionic metamaterials. Specifically, we perform high-throughput atomistic-spin-dynamic (ASD) simulations, using the UppASD package³⁸. The material class we propose consists of a monolayer film, which inherently supports non-trivial topologies^{17,39,40,41,42,43}. We present a theoretical pathway for constructing various complex topological textures from skyrmionic building-blocks, including reported high-order antiskyrmions and skyrmion bags, as well as the aforementioned skyrmionic metamaterials.

Results

Local assembly and interaction profiles

For the design of controlled assembly simulations in a 2D skyrmionic system, it is imperative to first analyze pairwise assembly, the relative stability of individual topological objects, and the interactions between the skyrmionic building blocks (i.e., between pairs of skyrmions, antiskyrmions, and skyrmioniums in all combinations). An efficient approach to study these behaviors involves examining the stationary configurations that emerge from performing spin dynamics simulations. There are two cases: some pairs attain a local energy minimum retaining a finite distance between them. For other pairs, the constituents will repel each other for large distances, but will merge (or annihilate) when forced closer to each other than a certain critical distance.

In our previous studies^17,44, we observed that Pd/Fe/Ir(111) can host a variety of high-order, small-sized skyrmionic textures, especially close to the transition region between skyrmion-lattice and ferromagnetic (single-domain) ground states. This agrees with several theoretical investigations that independently reported both antiskyrmions and higher-order topological spin configurations as metastable states in this system^{11,45,46,47,48,49}. In contrast to other systems that can host both skyrmions and antiskyrmions, such as Mn_1.4Pt_0.9Pd_0.1Sn⁵⁰, these small textures in Pd/Fe/Ir(111) are highly advantageous for the development of next-generation, high-information-density spintronics devices. Moreover, such sizes allow for a fully atomistic modeling with long-range magnetic interactions with ab-initio accuracy⁵¹. For these reasons, in the present paper, we focus our investigations on this system – but, as we discuss in subsequent sections, the ideas/results presented here may be extended to other frustrated materials; a detailed discussion of the generalization of results can be found in Section I E.

Unless otherwise mentioned, the Gilbert damping parameter in the ASD is α = 0.1 and we do all the simulations in a B^ext = 3.5 T magnetic field. We suggest this value of damping as a good lower-bound estimate for the α parameter of Pd/Fe/Ir(111), given the typical values for 3d ferromagnets/heavy-metal heterostructures⁵², and the enhancement influenced by the reduced thickness (single Fe monolayer)⁵³ and the presence of noncollinearity⁵⁴. For the chosen value of B^ext, the ground-state is ferromagnetic (single-domain).

To minimize thermal effects, not central to the study in this section, we set the simulation temperature to 1 mK. Generally, with fixed interaction parameters, a stronger applied field reduces both equilibrium and critical distances for combinations of skyrmions, antiskyrmions, and skyrmioniums, as shown in Fig. 1 a.

Design of 2D skyrmionic metamaterials through controlled assembly — **Fig. 1: Pair-wise interaction between topological building blocks.**

As shown in the top three subfigures of Fig. 1a, we observe distinct behaviors. The combination of skyrmion and skyrmionium (middle figure, top panel of Fig.1a) exhibits isotropic outer boundaries, resulting in rotation-invariant repulsion. However, the dynamics between the different systems presented in Fig. 1 differ. For instance, in skyrmion-skyrmion pairs, they separate until achieving an equilibrium distance through gyroscopic rotational dynamics⁵⁵.

In contrast, with skyrmion-skyrmionium pairs, only the skyrmionium moves while the skyrmion remains stationary. When an antiskyrmion, characterized by an anisotropic boundary, is part of the skyrmionic pair, the interaction behavior changes. Upon approaching the critical distance from far away, these components exhibit attraction, leading to assembly⁵⁶. As depicted in the bottom three subfigures of Fig. 1a, two skyrmionic textures assemble when their initial gap is smaller than the critical distance, resulting in either a high-order antiskyrmion or a vanishing ferromagnetic state. During the assembly of two antiskyrmions, as illustrated in Fig. 1b, a metastable “blocking point” emerges.

This object is in essence a quasi-equilibrium state that appears when two antiskyrmionic textures are close enough and the self-rotation is prevented in some way. Once a large enough lattice-like skyrmionic metamaterial is constructed, self-rotation becomes very difficult, since this would require the entire metamaterial to rotate. The result is that the lattice-like skyrmionic metamaterial becomes a very long-lived metastable state, despite its high complexity.

For a further discussion of the blocking-point concept, see Supplementary Note 2.

High-order antiskyrmion lattice

Three examples of self-stabilized lattice-like skyrmionic metamaterials that we have discovered are shown in the top subfigures of Fig. 2b–d). By using numerically prefabricated building blocks (see text below), i.e., high-order antiskyrmions with Q = − 2, − 5, and − 12, it is possible to construct such skyrmionic metamaterials, essentially similar to, e.g., artificial spin-ice structures and lattice based metamaterials. Here, high-order antiskyrmions represent magnetic textures with topological charge Q < − 1^17,57. To quantify the stability of the skyrmionic metamaterials of Fig. 2, we computed stability maps with temperature and external field as critical parameters. These stability maps are depicted at the bottom of Fig. 2. Note that in these figures stability is marked as a blue region (also marked with a 1) while the white region, marked with a 0, denotes an unstable configuration, that anihilates within the simulation time. In this stability analysis we terminated the simulations after 200 ps simulation time, and simply inspected if a configuration was stable by monitoring that it reached a static state. From this analysis, we can draw the conclusion that the magnetic configurations are long-lived (compared to, e.g., typical spin precession and relaxation timescales), with a lifetime exceeding 200 ps. Interestingly, the stability analysis of Fig. 2 shows that lattices of higher order skyrmions can survive significant temperature and field ranges. In this regard, it should be noted that a detailed analysis of the potential energy surface using geodesic nudged elastic band (GNEB)⁵⁸ could better quantify the lifetimes and transition energy barriers of the formed structures. However, applying GNEB to high-order topological structures in frustrated magnets remains challenging and will be explored in future work.

**Fig. 2: Construction and stability of lattice-like skyrmionic metamaterials.**

Below, we detail our approach for constructing and discovering long-lived lattice-like skyrmionic metamaterials. In essence, the first step is to construct high-order antiskyrmions (referred to as prefabricated building blocks henceforth) from simple antiskyrmions. Then, these prefabricated building blocks are combined together in such a way that self-rotation is prevented, which in turn stabilizes the resulting lattice-like skyrmionic metamaterial.

In this way we start by using single antiskyrmions to construct prefabricated building blocks via simulated controlled assembly. This prefabrication process is guided by the site-based topological charge distribution rule discussed in Section III and the fact that antiskyrmions show an attractive interaction with other structures (skyrmion, skyrmionium, antiskyrmion), as demonstrated in Section I A. A subtle issue in the prefabrication process is how to place the antiskyrmions so that they assemble to the building block we want to fabricate. In the simulation, this is done by placing them at the appropriate distances from each other.

In an experimental setting, this may be realized using, for example, gradient magnetic fields^{4,12,59,60,61}, electric fields created by the scanning tunneling microscope⁶², or even temperature gradients⁶³. In Fig. 2a we give one example on how to fabricate a Q = − 5 high-order antiskyrmion, with five controlled-assembly steps, starting from five isolated antiskyrmions (one at the center and four at the edges). During the simulation, we manually moved isolated antiskyrmions from the edges to the center and relaxed the spin structure to ensure the successful formation of the new high-order antiskyrmion and the completion of the fabrication process. Finally, the triangle-shaped high-order antiskyrmion (Q = −5) is fabricated by controlled placement of a single antiskyrmion. The energy change with each step is illustrated in the left part of Fig. 2a. This result indicates a dynamic construction method for building high-order antiskyrmions: sequentially manipulating pregenerated antiskyrmions. Once we have constructed the prefabricated building blocks, these can be used to construct lattice-like skyrmionic metamaterials. Numerous structures are of course possible, and which lattice-like skyrmionic metamaterial will result depends on the relative placement of the prefabricated building blocks. For instance, in order to create a lattice-like skyrmionic metamaterial with a hexagonal lattice as shown in Fig. 2b, the prefabricated building blocks are placed so that they form a hexagonal lattice.

Flake-like skyrmionic metamaterials

Together with the aforementioned lattice-like skyrmionic metamaterial, another possible stable type suggested in this work is the flake, which is shown in the middle subfigure of Fig. 3c and (d). To demonstrate the building pathway of a flake, we begin by constructing ring-like prefabricated building blocks (Q = −6). As shown on the left of Fig. 3a, according to the site-based topological charge distribution rule, there are two methods for forming this ‘ring’ structure. In the first method (top of Fig. 3a), six antiskyrmions (Q = −1) are positioned close enough, less than the critical distance, to the boundary of a skyrmionium (Q = 0) to ensure simultaneous attraction. A 300 ps spin dynamics simulation on this initial spin configuration predicts the assembly of the antiskyrmions, forming the target prefabricated building block with a total topological charge of −6. The second method (bottom of Fig. 3a) involves moving an antiskyrmion close to the notch of a crescent-shaped prefabricated building block (Q = −5), created from 5 antiskyrmions. In this case, only 100 ps ASD simulations are required to validate the formation of the target high-order antiskyrmion prefabricated building block.

After acquiring this high-order antiskyrmion, we use it as an elementary building block to construct a high-order antiskyrmion lattice, as shown in Fig. 3b. This skyrmionic lattice is predicted to be unstable even at 1 K, where a minor self-rotation can easily move the blocking point, leading to connections between components and unwinding. During this process, we also observe the formation of some flake-like textures.

Expanding on this, we increase the antiskyrmion density based on the previous building block. Ultimately, as shown in Fig. 3c, we isolate a honeycomb flake-like skyrmionic metamaterial (Q = −22). In Fig. 3d, by adding one more antiskyrmion, we obtain a honeycomb flake-like skyrmionic metamaterial (Q = −23) hosting the well-known 5-7 honeycomb lattice defects⁶⁴, where a 5-7 defect in a honeycomb lattice, is a defect in which one ring has 5 sites and an adjacent ring has 7 sites, instead of all rings having 6 sites. Both of those flake-like skyrmionic metamaterials are predicted to be stable even at 10 K in Pd/Fe/Ir(111). However, the flake with the 5–7 defect is less stable than the one without. To further analyze this, we define the relative defect formation energy as:

$${E}_{rdf}=frac{1}{n}left[{E}_{hf}+n({E}_{ASK}-{E}_{FM})-{E}_{df}right].$$

(1)

Here, E_rdf, E_hf, E_FM, E_ASK, and E_df represent the relative defect formation energy, the energy of the honeycomb flake-like skyrmionic metamaterial under a ferromagnetic background, the energy of the ferromagnetic background, the energy of an isolated antiskyrmion under a ferromagnetic background, and the energy of the honeycomb flake-like skyrmionic metamaterial with a 5-7 defect under a ferromagnetic background, respectively. In Eq. (1), n denotes the topological charge difference between the two flakes. Considering the DFT calculated magnetic spin moment of 2.78 μ_B/Fe atom, and the applied field B^ext = 3.5 T, the relative defect formation energy between the two flakes is E_rdf ~0.15 mRy/defect. This indicates that introducing defects to the flake-like skyrmionic metamaterial decreases its stability.

Stabilized cell-like skyrmionic metamaterial with texture nesting

To optimize the use of all basic building blocks relevant for this investigation, namely skyrmions, antiskyrmions, and skyrmioniums, as shown in Fig. 4e in this section we suggest several promising cell-like skyrmionic metamaterials and explore the construction pathway from small skyrmion bags to cell-like skyrmionic metamaterials, inspired by the design of skyrmion bags. These bags can also be interpreted as magnetic textures stabilized by skyrmionium. As demonstrated in Fig. 4a–d, we have been able to stabilize several skyrmions with different complexity and topological charges. This involves, e.g., skyrmion bags with topological charge ranging from −1 to −4, by applying a controlled assembly approach to sets of skyrmionium and antiskyrmions. Note that the information in Fig. 4 shows the time evolution of the different topological objects. At 0 ps we show the configuration obtained initially in the controlled assembly protocol, and the data at 200 ps (Fig. 4a–d) shows the final converged configuration. We observe that in the so conducted, controlled assembly process, the textures inherit an homogeneous, boundary from the skyrmionium, lose the property of attraction but retain their topological charge. Since these skyrmion bags also adhere to the local topological charge distribution, it is expected that more complex skyrmion bags can be constructed through similar operations.

**Fig. 4: Cell-like skyrmionic metamaterials with distinct topological charges.**

With the introduction of these new kinds of prefabricated building blocks, i.e., skyrmion bags, the repertoire of building blocks with isotropic boundary conditions, or in other words, those that exhibit mutual repulsion, is now expanded beyond the initial three types. Utilizing these building blocks, we introduce a texture nesting operation, which involves using antiskyrmions and skyrmionium to create line-like, closely-packed skyrmion bags as boundaries to nest other textures with isotropic boundary conditions. As shown in Fig. 4e, we conduct a 7 ns ASD simulation in a system of 300 × 300 spins. Due to the large system size and extended simulation time, we set the damping to 0.2 and the time step to 1 fs to expedite convergence. We design a texture, similar to but not identical to one in a previous study¹⁸, containing only skyrmions within the skyrmion-bag-ring. This results in a skyrmion bag with a skyrmion lattice inside. Subsequently, we experiment with embedding skyrmionium and small skyrmion bags within the larger skyrmion bag ring, as illustrated in the subsequent rows of Fig. 4e. This controlled assembly mechanism potentially leads to more complex metastable textures on demand. All these dynamic pathways can be summarized as utilizing repulsion from skyrmions and skyrmion bags to inhibit the shrinking process of a newly assembled skyrmion bag.

Generalization of results

Until now, all simulations have been performed using spin Hamiltonian parameters for Pd/Fe/Ir(111) calculated from Density Functional Theory (DFT)⁵¹. These calculations are known to give reliable values of, e.g., the Heisenberg exchange as well as the anisotropic Dzyaloshinskii-Moriya interaction (DMI). However, as with every theory, an uncertainty naturally exists – in part, due to the approximations involved in all theoretical steps and methods. In order to test the sensitivity of the results predicted here towards such uncertainties, we designed a set of simulations where the magnetic parameters are perturbed by given percentages. This also extends the controlled assembly protocols to encompass a wider range of materials. To achieve this and address more realistic conditions, we now intend to emulate a possibly common experimental situation: the assembly process occurring among spin textures surrounded by other metastable excitations^49,65,66. A simplified way to emulate this is to set up a relatively small (60 × 60) spin lattice with periodic boundaries (PBs). Hence, the set of simulations involves two skyrmioniums and one antiskyrmion in this lattice, and the application of a scalar perturbation to the DFT calculated parameters of Pd/Fe/Ir(111). Four obtained final states (1–4) after a 100 ps ASD simulation at each sampling point are shown in Fig. 5a subplots 1 to 4. In these simulations, as illustrated in Fig. 5b, we first applied random perturbations (± 8%) to the DFT-calculated DMI (D_ij), anisotropy (K^U), and Heisenberg exchange (J_ij) interactions.

**Fig. 5: Tests of the robustness and general applicability of the controlled assembly process.**

We found that the final state 4 is not strictly limited to the originally DFT calculated parameters.

Furthermore, to evaluate the general applicability of our findings, we choose to apply pairwise tuning of the coefficients in the spin Hamiltonian, e.g., the D_ij– J_ij pair. As shown in Fig. 5c, we first scaled the DFT calculated nearest-neighbor J_ij from 0.6 to 1.4 and varied the present damping (originally α = 0.1) from 0.1 to 0.9. The results indicate that with lower J_ij and higher damping, we can successfully form a skyrmion bag. However, increasing J_ij leads to more energetic spin wave emissions from the assembly process, which can destabilize structures like skyrmioniums and skyrmion bags compared to antiskyrmions. Interestingly, when J_ij is lower, 100 ps may not suffice for the skyrmion bag to transition from a sparse state (as shown in subplot 3 in Fig. 5a) to a dense state (as shown in subplot 4 in Fig. 5a). In turn, changes in the α parameter (see Fig. 5c) do not directly impact the final states, but instead influence on the lifetimes of spin-wave-like excitations and other transient modes emitted during the assembly process. Thus, for smaller α values these excitations can live longer, facilitating the system to escape shallower metastable states by transiently perturbing them; the probability that the system will stabilize in one of such metastable states decreases, but it is still possible, as shown in Fig. 5c for α = 0.2 or α = 0.3. This process is corroborated by the finite size of our simulation box with PBs, which eventually causes these excitations to reflect back to the higher-order skyrmionic textures and destabilize them (as in an area surrounded by other topological structures). With this insight, we reverted the Gilbert damping to its original setting (α = 0.1) to mitigate the effects of such spin-wave-like excitations. We then analyzed the dynamics of controlled assembly with perturbed interaction parameters (D_ij, J_ij, and K^U).

Similar generalization tests are performed in D_ij–K^U, J_ij–K^U, and D_ij–J_ij pairs, as shown in Figures 5d–f. From these results, it is clear that random perturbations in K^U do not significantly affect the formation of magnetic states. In turn, the behavior under perturbations in J_ij and D_ij is more complex: the final state depends, for example, on how the random perturbations affect the frustration in the Fe-Fe exchange interactions (the third-nearest shell in Pd/Fe/Ir(111), see Supplementary Fig. 1). Apart from this complex behavior, it is evident that the results are not limited to a fortuitous choice of the Pd/Fe/Ir(111) magnetic interaction parameters, but can survive random perturbations. Skyrmionic metamaterial and controlled assembly are thus anticipated to occur in other systems where long-range frustration and DMI interactions coexist.

Discussion

The ability to use simulation tools to pick up individual atoms and to place them where one wants, so that one can then investigate their possible interaction and stability, has been a fantastic asset in understanding chemical reactions, and therefore deepened the understanding of, e.g., catalysis, crystal growth and surface reconstructions. Several softwares have been designed for this task. The ability to do this with magnetic textures is far less investigated. Here we take steps to accommodate this, and we demonstrate that controlled assembly protocols, that combine artificial intelligence based segmentation tools with atomistic spin-dynamics simulations, offer a great opportunity for such investigations. In this work we began by systematically analyzing the pairwise interactions between skyrmionic elements, pinpointing their critical distances. This understanding facilitates the exploration of single high-order antiskyrmion formation through controlled assembly, guided by the distribution of site-based local topological charges. Although manipulating skyrmionic textures experimentally remains challenging, the controlled assembly method presents new potential solutions for constructing complex spin textures. Crucially, we successfully identified a variety of skyrmionic metamaterials. These structures are notable for their stability for a wide range of different Heisenberg and DM interactions, presenting new realms of complexity and stability in magnetic textures. The interactions used for the majority of the investigations presented here, represent a system that has been realized experimentally: the Pd/Fe/Ir(111) system. Hence, the spin Hamiltonian used here is appropriate for a system that has been observed in nature, as opposed to a model Hamiltonian where the interaction parameters may be tuned or fitted, having a specific finding or observation in mind. Judging from the results presented here, the Pd/Fe/Ir(111) system has a fantastic richness in it, in terms of the topological magnetic structures it can host.

Although the specific capabilities of identified skyrmionic metamaterials, as compared to conventional skyrmionic textures, require further elucidation, the findings of this study significantly contribute to the exploration of the theoretical landscape of skyrmionic textures. Specifically, our research reveals a new set of long-lived textures and their formation pathways, which could potentially act as novel candidate information carriers in the construction of spintronic devices. This work opens promising avenues for the practical realization of long-lived complex magnetic textures, i.e., constructing these textures from fundamental skyrmionic configurations. We anticipate that our work will catalyze further experimental explorations, thereby advancing the development of future spintronic devices.

Methods

We describe this system using a classical atomistic spin Hamiltonian, with all relevant interaction parameters derived from Density Functional Theory (DFT)^17,51:

$$begin{array}{lll}{mathscr{H}},=,-mathop{sum}limits_{ine j}{J}_{ij}{{bf{S}}}_{i}cdot {{bf{S}}}_{j}-mathop{sum}limits_{ine j}{{bf{D}}}_{ij}cdot left({{bf{S}}}_{i}times {{bf{S}}}_{j}right)\qquad;;; -mathop{sum}limits_{i}{mu }_{i}{{bf{B}}}^{{rm{ext}}}cdot {{bf{S}}}_{i}-mathop{sum}limits_{i}{K}_{i}^{{rm{U}}}{left({{bf{S}}}_{i}cdot {{bf{e}}}_{z}right)}^{2},end{array}$$

(2)

where S_i represents the normalized spin moment at site i, and μ_i its magnetic moment length. The parameters J_ij, D_ij, K^U, e_z, and B^ext correspond, respectively, to Heisenberg exchange interactions, Dzyaloshinskii-Moriya interactions, uniaxial anisotropy, the easy axis vector, and the applied magnetic field. The DFT calculated parameters are scaled and perturbed in the generalization section using preset factors and restricted Gaussian noise (see Supplementary Note 1). This is done to enhance the generalizability and robustness of the model. The stability test, involving perturbation of input parameters, demonstrates that our findings are not limited to the Pd/Fe/Ir system but also applicable to a broader range of skyrmion-hosting thin-film materials. While dipole-dipole interactions (DDI) are known to stabilize topological skyrmionic structures^50,66,67, including an explicit DDI term in the spin Hamiltonian (Eq. (2)) had minimal impact on the results (see discussion in Supplementary Note 3).

Using the spin Hamiltonian of Eq. (2), we performed high-throughput ASD simulations in a wide range of initial spin configurations. The ASD simulations were carried out in for a single fcc(111) monolayer (i.e., a hexagonal lattice) representing the Fe layer in the ideal fcc-Pd/Fe/Ir(111) structure, with the Ir experimental lattice parameter (Fe-Fe equilibrium distance of 0.272 nm). The initial spin configurations, where the skyrmionic textures are separated by various distances, can be used to understand their interactions. To facilitate the construction of systems involving topological structures separated by controlled distances, as shown in Fig. 6a, we employed a state-of-the-art artificial intelligence segmentation tool to extract elemental spin textures from the optimized larger spin system^17,68. This tool allows for the systematic design of initial spin textures, such as selecting different element textures, altering the distances between elements, and rotating element textures. We studied the motion of skyrmionic textures at T = 1 mK using atomistic spin dynamics, numerically solving the stochastic Landau-Lifshitz-Gilbert equation as implemented in the Uppsala Atomistic Spin Dynamics (UppASD) package^38,69:

$$frac{d{{bf{S}}}_{i}}{dt}=-{gamma }_{{rm{L}}}{{bf{S}}}_{i}times {{bf{B}}}_{i}^{{rm{eff}}}(t)-{gamma }_{{rm{L}}}alpha {{bf{S}}}_{i}times left({{bf{S}}}_{i}times {{bf{B}}}_{i}^{{rm{eff}}}(t)right).$$

(3)

**Fig. 6: Design of initial spin configurations and site-based local topological charge density.**

In Eq. (3), α is the scalar Gilbert damping constant, ({gamma }_{L}=frac{gamma }{left(1+{alpha }^{2}right)}) is the renormalized gyromagnetic ratio, and ({{bf{B}}}_{i}^{{rm{eff}}}(t)) is the effective magnetic field at the site i:

$${{bf{B}}}_{i}^{{rm{eff}}}(t)={{bf{B}}}_{i}+{{bf{B}}}_{i}^{{rm{f}}}(t)=-frac{1}{{mu }_{i}}frac{partial {mathscr{H}}}{partial {{bf{S}}}_{i}}+{{bf{B}}}_{i}^{{rm{f}}}(t),$$

(4)

which here includes a time-dependent stochastic thermal field ({{bf{B}}}_{i}^{{rm{f}}}(t)). The random field ({{bf{B}}}_{i}^{{rm{f}}}) is assumed to be a Gaussian stochastic process that satisfies the time-averages (langle {B}_{i}^{{rm{f}},alpha }(t){B}_{j}^{{rm{f}},beta }({t}^{{prime} })rangle =2{D}_{i}{delta }_{ij}{delta }_{alpha beta }delta (t-{t}^{{prime} })) and (langle {B}_{i}^{{rm{f}},alpha }rangle =0), where α, β are Cartesian coordinates, δ_ij is the Kronecker delta, (delta (t-{t}^{{prime} })) is the Dirac delta, and D_i = αk_BT/(γμ_i) is the strength of the thermal fluctuations, determined by the fluctuation-dissipation theorem³⁸.

As depicted in Fig. 6b and c, we propose a protocol to guide the construction of all skyrmionic metamaterials in 2D films, utilizing a site-based topological charge distribution method. In this context, a ’site’ is a qualitative idea that refers to a representative piece of the texture, characterized by a unique local topological charge density distribution. A total of four distinct sites have been defined in this work, with more detailed information available in the Supplementary Note 2.

The topological charge Q of a given continuous two-dimensional magnetization texture with unitary local direction (overrightarrow{n}(x,y)) is computed using

$$Q=frac{1}{4pi }{int}_{A}overrightarrow{n}cdot left(frac{partial overrightarrow{n}}{partial x}times frac{partial overrightarrow{n}}{partial y}right)dxdy,$$

(5)

where the integral is taken over the area A in the xy-plane. Similarly, for a discrete spin system, one can use the Berg and Lüscher invariant^70,71 which reads

$$Q=sum _{{x}^{* }}q({x}^{* }),$$

(6)

for all unique topological densities q(x^*) attached to vertices x^* in the lattice with spins sitting around, also taking into account the periodic boundary conditions. By construction, both Eqs. (5) and (6) result in the total topological charge existent in the system. The discrete format allows to construct the concept of a local topological charge, q_loc, which can be defined as the following summation:

$${q}_{{rm{loc}}}=sum _{{x}^{* }in {x| r(x)le {r}_{c}}}q({x}^{* }),$$

(7)

where r(x) denotes the distance of the point x from a reference point in the lattice, and r_c represents a cutoff radius⁷¹. In this case, obviously (mathop{lim }nolimits_{{r}_{c}to infty }{q}_{{rm{loc}}}=Q).