[1]Kyoung-Min, Kim [2]Se Kwon, Kim

[1] Center for Theoretical Physics of Complex Systems, Institute for Basic
Science, 55 Expo-ro, Yuseong-gu, 34126, Daejeon, Republic of Korea

[2] Department of Physics, Korea Advanced Institute of Science and Technology,
291 Daehak-ro, Yuseong-gu, 34141, Daejeon, Republic of Korea

3] Department of Physics, Hanyang Univercity, 222 Wangsimni-ro, Seongdong-gu,
04763, Seoul, Republic of Korea

# Emergence of stable meron quartets in twisted magnets

<EMAIL_ADDRESS>Gyungchoon  Go  Moon Jip  Park<EMAIL_ADDRESS>*  *  [

###### Abstract

The investigation of twist engineering in easy-axis magnetic systems has
revealed the remarkable potential for generating topological spin textures,
such as magnetic skyrmions. Here, by implementing twist engineering in easy-
plane magnets, we introduce a novel approach to achieve fractional topological
spin textures such as merons. Through atomistic spin simulations on twisted
bilayer magnets, we demonstrate the formation of a stable double meron pair in
two magnetic layers, which we refer to as the “Meron Quartet” (MQ). Unlike
merons in a single pair, which is unstable against pair annihilation, the
merons within the MQ exhibit exceptional stability against pair annihilation
due to the protective localization mechanism induced by the twist that
prevents the collision of the meron cores. Furthermore, we showcase that the
stability of the MQ can be enhanced by adjusting the twist angle, resulting in
increased resistance to external perturbations such as external magnetic
fields. Our findings highlight the twisted magnet as a promising platform for
investigating the intriguing properties of merons, enabling their realization
as stable magnetic quasiparticles in van der Waals magnets.

###### keywords:

van der Waals magnet, moiré magnet, twist engineering, magnetic vortex, meron,
topological spin texture

###### Contents

  1. 1 Introduction
  2. 2 Results
    1. 2.1 Antiferromagnetic domain array
    2. 2.2 Emergence of stable merons: meron quartets
    3. 2.3 Stability of meron quartet states
    4. 2.4 Diverse forms of meron quartets
  3. 3 Discussion
  4. 4 Methods
    1. 4.1 Moiré superlattice
    2. 4.2 Interlayer Heisenberg exchange interactions
    3. 4.3 Iterative optimization method
    4. 4.4 Local magnetic energy maps
    5. 4.5 Determination of FM-MD phase boundary
    6. 4.6 Critical twist angle formula for FM-MD transition
    7. 4.7 Determination of MD-MQ phase boundary
    8. 4.8 Determination of critical field strengths for MQ-MD transition
    9. 4.9 MQ states with different skyrmion numbers
  5. 5 Supplementary information
    1. 5.1 Supplementary Video 1
    2. 5.2 Supplementary Video 2

###### List of Figures

  1. 1 Schematic illustration
  2. 2 Emergence of AFM domain array
  3. 3 Emergence of stable meron-antimeron pairs
  4. 4 Comparison of MQ and MD
  5. 5 Stability of MQ
  6. 6 Ext. Data: Interlayer exchange interactions
  7. 7 Ext. Data: Relaxation method
  8. 8 Ext. Data: FM-MD transition
  9. 9 Ext. Data: MD-MQ transition
  10. 10 Ext. Data: Critical field for MQ-MD transition
  11. 11 Ext. Data: MQ states with different skyrmion numbers

##  1 Introduction

Figure 1: Schematic illustration of a stable meron-antimeron pair in a twisted
magnet. a: Twist-induced antiferromagnetic (AFM) domain array in a
ferromagnetic (FM) order background. Red and blue colors indicate parallel and
antiparallel spin alignments between the top and bottom layers, respectively.
b: Emergence of a stable meron-antimeron pair. Arrows and circles depict their
in-plane winding textures and core positions, respectively. c: Schematic
energy landscape for localizing the cores within different AFM domains (left),
and its schematic profile along the dashed line (right).

The search for novel spin configurations motivated by fundamental interest and
technological applications has led to the discovery of intriguing textures
with nontrivial topology in various magnetic systems 1. Specifically, vortex-
type topological spin textures, so-called merons, have been observed in
confined magnetic disks 2, 3, 4, 5, 6 and continuous thin films 7, 8, 9, 10,
11, 12, 13. Recently, monolayer chromium trichloride (CrCl3), a two-
dimensional (2D) van der Waals (vdW) magnetic crystal, has emerged as a
promising candidate for achieving merons in this novel atomically thin limit
14, 15, 16. The intrinsic easy-plane magnetic anisotropy in such a system
offers a pathway to attain in-plane swirling spin textures for merons.
However, these merons are inherently unstable against pair annihilation and
possess only a limited lifespan 16. The mechanism responsible for stabilizing
such merons has yet to be definitively determined, hampering future
exploration of merons in 2D vdW magnets.

The field of twist engineering has opened up a fascinating realm of
possibilities in generating topological spin textures in 2D vdW magnets. By
harnessing moiré patterns, researchers have demonstrated the creation of
skyrmion spin textures in lattice-mismatched heterostructures 17, 18, 19 and
twisted homo-bilayer systems 20, 21, 22, 23, 24, 25, with a primary focus on
Ising-type easy-axis magnetic systems. However, extending this approach to XY-
type easy-plane magnets holds tremendous intrigue, as it enables exploration
of captivating phenomena such as the Berezinskii-Kosterlitz-Thouless
transition 26, 27, the emergence of merons 15, 16, and the potential discovery
of hidden magnetic phases driven by strong spin fluctuations 28. The recent
discovery of twisted magnets further enhances the allure, prompting a
continued investigation into the captivating moiré effects in these systems
29, 30, 31, 32.

In this study, we investigate twist engineering in vdW magnets as a promising
avenue to realize stable merons. By conducting atomistic spin simulations, we
demonstrate that antiferromagnetic domain arrays in the twisted magnet 33, 23,
24, 29, 30, 31, 32, 34 can be utilized to localize the cores of merons along
the boundaries of their respective domains (Fig. 1a-b). This localization
mechanism effectively preserves the stable spin configuration of the meron
pair by separating their cores (Fig. 1c). Furthermore, we show that the
stability of merons can be tuned by adjusting the twist angle, providing
controllable resistance against external magnetic fields. These findings
present a promising avenue for achieving stable merons as magnetic
quasiparticles in vdW magnets, offering the opportunity to explore their
captivating properties with significant flexibility through external stimuli
35, 36 or the creation of heterostructures 37.

##  2 Results

###  2.1 Antiferromagnetic domain array

Figure 2: Emergence of AFM domain array. a: Moiré superlattice for a twist
angle $\theta=5.08$°. The colored circles denote local stacking patterns,
including AA (green), AB (blue), BA (cyan), and monoclinic (red). The yellow
rhombus and black arrows denote the unit cell and lattice vectors of the moiré
superlattice, respectively. b: Modulation of the local interlayer exchange
energy ($J_{i}^{\perp}$) computed for $\theta=1.02$°, with blue (red) color
representing FM (AFM) coupling. c: Zero temperature magnetic phase diagram
depicting FM phase (white) and magnetic domain phase (green), as a function of
twist angle ($\theta$) and intralayer exchange ($J$). The markers represent
the phase boundary obtained from numerical simulations, while the dotted line
indicates the fitting curve computed from an effective continuum model (see
Methods 4.6). d-f: Ground-state spin configuration in the magnetic domain
phase for $J=2$ meV and $\theta=1.61$°. d-e: The color scales denote the phase
angles ($\phi_{t,b}$) of the normalized spin vectors
$\mathbf{n}_{t,b}=(\cos{\phi_{t,b}},\sin{\phi_{t,b}},0)$ in the top (d) and
the bottom (e) layers, respectively. In the magnified images, the arrows
denote the direction of $\mathbf{n}_{t}$ and $\mathbf{n}_{b}$, and the yellow
areas highlight the domain walls with $90$° or $-90$° spin rotations. Here,
“t” and “b” represent the top and bottom layers, respectively. f: The color
scale denotes the relative orientation of the spin vectors between the two
layers ($\mathbf{n}_{t}\cdot\mathbf{n}_{b}$), where red (blue) represents
parallel (antiparallel) alignment. The magnified image highlights the
correspondence between the interlayer exchange and magnetic domain structure.

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

$\displaystyle H=$ $\displaystyle-\frac{J}{2}\sum_{l=t,b}\sum_{\langle
i,j\rangle}\mathbf{S}_{i}^{l}\cdot\mathbf{S}_{j}^{l}+A\sum_{l=t,b}\sum_{i}\big{(}\bm{S}_{i}^{l}\cdot\hat{z}\big{)}^{2}$
$\displaystyle+\sum_{i,j}J_{ij}^{\perp}\mathbf{S}_{i}^{t}\cdot\mathbf{S}_{j}^{b}.$
(1)

Here, $\mathbf{S}_{i}^{l}$ represents the spin at site $i$ on the top layer
($l=t$) and the bottom layer ($l=b$). $J$ represents the intralayer FM
exchange interactions between nearest-neighbor spins. $A=0.1$ meV represents
the single-ion anisotropy energy favoring in-plane magnetization.
$J_{ij}^{\perp}$ represent the interlayer exchange interactions, which are
adopted from previous ab-initio calculations on bilayer CrI3 24. Due to so-
called stacking-dependent interlayer magnetism 38, 39, 40, 41, 42, 21,
$J_{ij}^{\perp}$ switch from FM to AFM coupling depending on the local
stacking pattern between the two magnetic layers (Methods 4.2). Consequently,
the interlayer exchange coupling exhibits the coexistence of AFM and FM
interactions in the moiré superlattice accommodating various local stacking
patterns (Fig. 2a)20, 21, 22, 23, 24, 25, 34. We illustrate this behavior in
Fig. 2b through the map of the local interlayer exchange energy
$J_{i}^{\perp}=\sum_{j}J_{ij}^{\perp}$ computed in an FM configuration
$\mathbf{S}_{i}^{l}=S\hat{z}$ 20, 22, 24, 34. Specifically, in the monoclinic
stacking region (red patches), $J_{i}^{\perp}$ exhibits AFM character
($J_{i}^{\perp}>0$), indicating a tendency for the spins in the top and bottom
layers to align antiparallel to each other. Conversely, in the other stacking
regions, $J_{i}^{\perp}$ exhibits FM character ($J_{i}^{\perp}<0$), signifying
a preference for parallel alignment. As a result, the twisted magnet embeds
local AFM patches in a background of FM coupling (Fig. 2b).

In this work, we investigate the influence of the AFM patches on the
stabilization of merons. We first identify the magnetic phases of twisted
easy-plane magnets. Our atomistic simulations on Eq. (2.1), conducted using an
iterative optimization method (Methods 4.3), reveal the zero-temperature
magnetic phase diagram shown in Fig. 2c. Within this diagram, we observe two
distinct magnetic phases: an FM phase and a magnetic domain (MD) phase. The FM
phase exhibits a uniform spin configuration with parallel alignment between
the spins of the top and bottom layers. On the other hand, the MD phase
exhibits antiparallel alignment within the AFM patches, while maintaining
parallel alignment outside these patches to minimize the interlayer exchange
energy (Fig. 2d-f). This contrasting AFM-FM order results in the formation of
AFM domains within each AFM patch as well as domain walls surrounding the
domains, characterized by spin rotations of $90$° and $-90$°. Furthermore, the
AFM domains are arranged into an array structure across the superlattice,
resembling a Kagome lattice. We dub this distinctive magnetic structure an AFM
domain array.

We attribute the emergence of the AFM domain array to the amplified effect of
interlayer exchange in the small twist angle regime
($\theta<\theta_{c1}\sim\sqrt{\bar{J}_{\perp}/J}$). In this regime, the
formation of an AFM domain becomes energetically favorable as the reduction in
the interlayer exchange energy ($\Delta
E_{\textrm{inter}}\sim-\bar{J}_{\perp}\frac{L^{2}}{a^{2}}\sim-\bar{J}_{\perp}\theta^{-2}$)
outweighs the increase in the intralayer exchange energy ($\Delta
E_{\textrm{intra}}\sim J$). Here, we consider $\bar{J}_{\perp}$ as the average
value of $J^{\perp}_{i}$ within the AFM patch. The patch is approximated as a
disk with a radius denoted by $L\sim a/\theta$, where $a$ represents the
lattice constant of the honeycomb lattice. Despite interlayer exchange being
weaker than intralayer exchange, its effect is significantly amplified by the
large size of the AFM patch ($\pi L^{2}$). Consequently, the formation of the
AFM domain array is expected when the twist angle is sufficiently small. This
phenomenon manifests in the phase diagram, which exhibits a consistent
relationship $\theta_{c1}\sim\sqrt{\bar{J}_{\perp}/J}$ for the phase boundary
between the FM and MD phases (dashed line in Fig. 2c).

###  2.2 Emergence of stable merons: meron quartets

Figure 3: Snapshots depicting the relaxation process and the formation of
stable meron-antimeron pairs. a-b: Intermediate state illustrating the
spontaneous generation of merons ($\textrm{M}_{\textrm{t1}}$,
$\textrm{M}_{\textrm{t2}}$, $\textrm{M}_{\textrm{b1}}$, and
$\textrm{M}_{\textrm{b2}}$) and antimerons
($\overline{\textrm{M}}_{\textrm{t1}}$, $\overline{\textrm{M}}_{\textrm{t2}}$,
$\overline{\textrm{M}}_{\textrm{b1}}$, and
$\overline{\textrm{M}}_{\textrm{b2}}$). c-d: Fully-relaxed state displaying
stabilized meron-antimeron pairs
($\textrm{M}_{\textrm{t1}}$–$\textrm{M}_{\textrm{t1}}$,
$\textrm{M}_{\textrm{b1}}$–$\textrm{M}_{\textrm{b1}}$), with the annihilation
of pairs $\textrm{M}_{\textrm{t2}}$–$\overline{\textrm{M}}_{\textrm{t2}}$ and
$\textrm{M}_{\textrm{b2}}$–$\overline{\textrm{M}}_{\textrm{b2}}$. Panels a,c
and b,d represent the top and bottom layers, respectively, corresponding to a
magnified area shown in Fig. 2d-e. Arrows indicate in-plane components, while
the color scale in the markers represents out-of-plane components. Marker
sizes are adjusted for better visibility. Shaded areas indicate AFM patches.
The parameters $J=2$ meV and $\theta=1.61$° are utilized.

A meron is a vortex-like topological spin texture, characterized by an integer
winding number. In contrast to conventional magnetic vortices, the core of a
meron exhibits an out-of-plane polarization, resulting in a unique half-
skyrmion number denoted as 12:

$\displaystyle Q=\frac{1}{4\pi}\int dx\int
dy(\partial_{x}\mathbf{n}\times\partial_{y}\mathbf{n})\cdot\mathbf{n}=p\cdot
w=-\frac{1}{2}.$ (2)

Here, the vector field
$\mathbf{n}=(\sin{\vartheta}\cos{\varphi},\sin{\vartheta}\sin{\varphi},\cos{\vartheta})$
represents the orientation of spins. The polarity $p$ and vorticity $w$ of
$\mathbf{n}$ are defined by
$p=\frac{1}{2}[\cos\vartheta(r=\infty)-\cos\vartheta(r=0)]$ and
$w=\frac{1}{2\pi}\oint_{\gamma}d\mathbf{l}\cdot\nabla\varphi$, respectively,
where $\gamma$ is any contour that encircles the core. Merons possess two
distinct characteristics, corresponding to the combinations
$(w,p)=\left\\{(+1,-\frac{1}{2}),(-1,+\frac{1}{2})\right\\}$. Antimerons are
counterparts to merons, possessing an opposing skyrmion number of
$Q=+\frac{1}{2}$ with two distinct characteristics
$(w,p)=\left\\{(-1,-\frac{1}{2}),(+1,+\frac{1}{2})\right\\}$.

In continuous magnetic systems, merons and antimerons typically exist as pairs
with opposite winding numbers ($w=+1$ and $w=-1$) 7, 8, 9, 10, 11, 12, 13. The
formation of such pairs allows their swirling spin textures to cancel out away
from the cores, resulting in localized spin configurations with finite energy.
However, the mutual attraction between the cores renders the pairs inherently
unstable, leading to pair annihilation during magnetization dynamics 43.
Consequently, in conventional untwisted magnetic systems, these magnetic
textures are usually observed as transient states with a limited lifespan 7,
8, 13, 16.

In this work, we discover that the merons and antimerons in twisted magnets
can evade pair annihilation by forming a double meron pair in two magnetic
layers. To illustrate the emergence of such stable merons, we present a
typical relaxation process of the magnetic state in Fig. 3, which is obtained
through the relaxation of a random initial configuration (Methods 4.7). In the
intermediate magnetic state (Fig. 3a-b), we observe the spontaneous formation
of four merons ($\textrm{M}_{\textrm{t1}}$, $\textrm{M}_{\textrm{t2}}$,
$\textrm{M}_{\textrm{b1}}$, and $\textrm{M}_{\textrm{b2}}$) and four
antimerons ($\overline{\textrm{M}}_{\textrm{t1}}$,
$\overline{\textrm{M}}_{\textrm{t2}}$, $\overline{\textrm{M}}_{\textrm{b1}}$,
and $\overline{\textrm{M}}_{\textrm{b2}}$) on the top and bottom layers. Upon
subsequent relaxation, the intra-patch pairs, i.e. the meron-antimeron pairs
occupying the same AFM patch within the same layer such as
$\textrm{M}_{\textrm{t2}}$–$\overline{\textrm{M}}_{\textrm{t2}}$ and
$\textrm{M}_{\textrm{b2}}$–$\overline{\textrm{M}}_{\textrm{b2}}$, undergo pair
annihilation due to the attractive interactions driven by the intralayer
exchange interactions. However, the inter-patch pairs, i.e. the meron-
antimeron pairs occupying different AFM patches within the same layer such as
$\textrm{M}_{\textrm{t1}}$–$\overline{\textrm{M}}_{\textrm{t1}}$ and
$\textrm{M}_{\textrm{b1}}$–$\overline{\textrm{M}}_{\textrm{b1}}$, remain
robust against pair annihilation due to their mutual correlation facilitated
by the interlayer coupling. As a result, the fully relaxed state (Fig. 3c-d)
accommodates only this correlated double meron pair.

We attribute the stabilization of the double meron pair to the localization of
their cores. During the relaxation process (Supplementary Video 1), we observe
that the cores shift exclusively along the boundaries of the AFM patches. This
behavior arises from the bulk energy minimization condition: the requirement
to minimize the interlayer exchange energy over the bulk region by maintaining
the AFM domain configuration, i.e. antiparallel and parallel alignments inside
and outside the AFM patches as depicted in Fig. 2f. Enforcing such a condition
localizes the cores within their respective AFM patches by constraining their
motions along the patch boundaries and prohibiting them from transferring to
other patches. This localization mechanism preserves the inter-patch pairs by
ensuring the separation of their cores and protecting them against pair
annihilation. However, the localization mechanism cannot preserve the intra-
patch pairs, as their cores reside within the same patch (Supplementary Video
2).

The localization mechanism of the meron pair can be understood by considering
the effective confining potential arising from the interlayer coupling (Fig.
1c). The cores of merons in one layer (e.g., the top layer) experience an
effective potential generated by the other layer (e.g., the bottom layer)
through the interlayer exchange coupling (Fig. 4m). This potential reaches its
minimum energy along the boundaries of the AFM patches due to the bulk energy
minimization condition. Any displacement of the cores away from these
boundaries results in increased energy relative to the minimum, creating
potential wells along the AFM patch boundaries. These wells act as confining
forces, effectively localizing the cores within them. Furthermore, the
establishment of such potential wells requires the presence of two counterpart
merons in the bottom layer to facilitate the bulk energy minimization
condition. Consequently, the creation of four merons is required to protect
the merons via the confining potential.

Figure 4: Comparison of meron quartet (MQ) state to magnetic domain (MD)
state. a-c/d-f: Spin configurations of the MQ (a-c) and MD (d-f) state,
corresponding to Fig. 3b,d and Fig. 2d-f, respectively. a-b/d-e: Phase angles
($\phi_{t,b}$) of the normalized spin vectors
($\mathbf{n}_{t,b}=(\cos{\phi_{t,b}},\sin{\phi_{t,b}},0)$) in the top (a/d)
and bottom (b/e) layers. In a-b, the black color indicates out-of-plane
polarization. c/f: Relative orientation ($\mathbf{n}_{t}\cdot\mathbf{n}_{b}$)
between the top and bottom layers, with red (blue) denoting parallel
(antiparallel) alignment. g-i/j-l: Local magnetic energy maps corresponding to
the spin configurations (a-c/d-f). g-h/j-k: Intralayer exchange energy plus
single-ion anisotropy energy ($\textrm{E}_{\textrm{t,b}}$) in the top (g/j)
and bottom (h/k) layers, respectively. i/l: Interlayer exchange energy
($\textrm{E}_{\textrm{tb}}$). In a-i, the dashed lines denote the boundaries
of AFM patches. m/n: Schematic illustration of the MQ (m) and MD (n) states
across a single AFM patch, corresponding to the dotted lines shown in a-b and
d-e, respectively. Red and blue colors depict the AFM patch region
($J_{i}^{\perp}>0$) and FM coupling background ($J_{i}^{\perp}<0$),
respectively. Yellow and blue arrows represent spin orientations in the top
and bottom layers, respectively.

Based on our findings, we introduce a novel magnetic state dubbed the “Meron
Quartet” (MQ) state, which consists of four merons, two for each layer, as
depicted in Fig. 4. This state exhibits two key characteristics: Firstly, each
layer contains two vortices with opposite winding numbers ($w=+1$ and $w=-1$),
ensuring the total winding number cancels out. Secondly, each occupied AFM
patch harbors two vortices with the same winding numbers ($w=+1$ or $w=-1$) in
both the top and bottom layers. These specific arrangements enable the MQ
state to realize stable meron pairs through the implementation of the
confining potential. Furthermore, the MQ state minimizes the interlayer
exchange energy over the bulk region (Fig. 4c,i), as observed in the ground MD
state (Fig. 4f,l). The distinction between the MQ and MD states lies solely in
the localized core energy (Fig. 4g-h vs. j-k). As a result, the MQ state
exhibits a low magnetic energy of -3.064 meV per spin, comparable to the
energy of the ground MD state (-3.071 meV), despite its intricate spin
textures arising from the presence of merons (Fig. 4a-b). In other words, the
MQ state naturally accommodates the stable correlated pairs of merons, i.e.
the meron quartet, while satisfying the bulk energy minimization condition
(Fig. 4m), similar to the ground MD state (Fig. 4n).

###  2.3 Stability of meron quartet states

Figure 5: a: Magnetic phase diagram illustrating the meron quartet phase
(green) and the magnetic domain phase (white), as a function of twist angle
($\theta$) and intralayer exchange ($J$). The markers represent the phase
boundary determined through numerical simulations, while the dashed line
represents the phenomenological fitting curve $\theta_{c2}=\sqrt{J_{c2}/J}$,
where the fitting parameter $J_{c2}$ is found to be $8$ meV. b: Out-of-plane
net magnetization ($M_{z}=\frac{1}{N}\sum_{l=t,b}\sum_{i}n_{i,z}^{l}$) as a
function of an external magnetic field in the out-of-plane direction ($B$),
with the solid and dashed lines corresponding to the MQ and MD states, as
depicted in Fig. 4a-c and d-f, respectively. The arrows mark the critical
field strengths ($B_{c1}$ and $B_{c2}$) that signify the degradation of the MQ
state to the MD state, as shown by the merging of the two curves at $B_{c1}$
and $-B_{c2}$. c: Evolution of $B_{c1}$ and $B_{c2}$ with twist angle
($\theta$). Error bars denote standard errors from different samples. The
parameter $J=2$meV is utilized.

We observe that the MQ state acquires metastability in the small twist angle
regime (green area in Fig. 5a). This phenomenon is attributed to the
enhancement of the confining potential in such a regime. To elucidate this, we
consider the transition of the MQ state to the MD state. This transition
necessitates the mutual attractions of the meron cores. However, such
attraction leads to an inevitable increase in the interlayer exchange energy,
as it disrupts the bulk energy minimization condition, as illustrated in Fig.
4c,i. We estimate this energy increase as $\Delta E_{\textrm{tb}}\sim
J_{\perp}^{\textrm{FM}}\frac{dL}{a^{2}}$, where $J_{\perp}^{\textrm{FM}}$
represents the FM interlayer exchange near an AFM patch, and $d$ signifies the
displacement of a core from its equilibrium position due to the attraction.
The term $\frac{dL}{a^{2}}\sim\frac{d}{a\theta}$ corresponds to the number of
spins that undergo unfavorable ordering for the interlayer exchange, which
increases as the twist angle $\theta$ decreases. Consequently, this energy
cost escalates as the twist angle $\theta$ decreases, creating a substantial
energy barrier for the transition to the MD state, which corresponds to the
confining potential mentioned before.

This high energy barrier surpasses the attractive force between merons and
antimerons, stabilizing the MQ state. The attractive force is mediated by the
intralayer exchange interaction energy between a meron and an antimeron within
the same layer. This energy is akin to the Coulomb energy between a vortex and
an antivortex in the XY model, denoted as $E_{\textrm{C}}\sim J\ln{(R/a)}$ 43,
where $R\sim a/\theta$ represents the distance between the cores. The
attraction of the cores by a displacement of $d$ can potentially reduce the
Coulomb energy by $\Delta E_{\textrm{C}}\sim-Jd/R\sim-Jd\theta/a$. However,
this energy reduction $\Delta E_{\textrm{C}}$ is surpassed by the energy
increase $\Delta E_{\textrm{tb}}$ in the small twist angle regime
($\theta<\theta_{c2}\sim\sqrt{J_{\perp}^{\textrm{FM}}/J}$). As a result, the
attraction of the cores is effectively prohibited in such a regime, leading to
the stabilization of the MQ state. We find that this phenomenon is evident in
the phase diagram, which exhibits a consistent relationship
$\theta_{c2}=\sqrt{J_{c2}/J}$ (dashed line in Fig. 5a), where $J_{c2}$ is a
fitting parameter proportional to $J_{\perp}^{\textrm{FM}}$.

The high energy barrier also indicates the enhanced stability of the MQ state
against external perturbations, such as external magnetic fields, in the small
twist angle regime. To further illustrate this stability, we incorporate the
Zeeman term:

$\displaystyle H_{\textrm{Zeeman}}=-g\mu_{B}B\sum_{l=t,b}\sum_{i}S_{i,z}^{l},$
(3)

where $B$ represents an external magnetic field applied in the out-of-plane
direction. Through the systematic examination of the behavior of the MQ state
under the influence of the Zeeman term (Methods 4.8), we identify the critical
field strength for the destruction of the MQ state, as illustrated in Fig. 5b.
The critical field strength differs depending on the relative orientation of
the applied field with respect to the polarity of the meron and antimeron
cores, with the antiparallel field exhibiting a much higher critical field
strength ($B_{c2}$) compared to the parallel field ($B_{c1}$). Nevertheless,
both critical field strengths are significantly enhanced as the twist angle
decreases (Fig. 5c). This corroborates the enhanced stability of the MQ state
in the small twist angle regime.

###  2.4 Diverse forms of meron quartets

We discover that the twisted magnet can realize stable merons in diverse forms
with a different total skyrmion number ($Q_{\textrm{tot}}=0,\pm 1,\pm 2$) and
distinct combinations of merons and antimerons between two magnetic layers in
addition to the specific configuration shown in Fig. 4a-c. Due to the
condition of the vanishing of the total winding number, two vortices
consisting of each inter-patch pair must exhibit opposite winding numbers
$w=+1$ and $w=-1$, respectively. However, the cores have the flexibility to
possess their own polarity, which can be either $p=-\frac{1}{2}$ or
$p=+\frac{1}{2}$. This gives rise to four potential configurations for each
pair: (i) M–$\overline{\textrm{M}}$, (ii) $\overline{\textrm{M}}$–M, (iii)
M–M, and (iv) $\overline{\textrm{M}}$–$\overline{\textrm{M}}$. Furthermore,
these configurations can independently occur in the top and bottom layers,
resulting in a total of sixteen potential configurations for the MQ state,
with distinct numbers of merons and antimerons in each layer. Our
investigation employing general random initial configurations confirms the
emergence of such diverse configurations (Fig. 10). This observation
highlights the flexibility in achieving merons and antimerons in the twisted
magnet, setting it apart from conventional magnetic systems, which typically
have fixed skyrmion numbers 1, 13.

##  3 Discussion

We have shown that the AFM domain array induced by the twist provides a
favorable environment for hosting stable merons. Once these topological
defects form, their destruction is impeded due to the constraints imposed by
the bulk energy minimization condition. Moreover, our theory suggests the
feasibility of realizing such stable merons and their enhanced robustness in
the small twist angle regime.

We propose our theory can be applied to CrCl3, which exhibits two essential
factors for hosting the meron quartet: easy-plane magnetic anisotropy 15 and
stacking-dependent interlayer magnetism 21. Notably, merons have been observed
in monolayer CrCl3 15, 16. In this context, twist engineering offers an
effective approach to realizing the meron quartet. Another potential
application lies in CrI3, where the modification of its intrinsic easy-axis
magnetic anisotropy to easy-plane anisotropy can be achieved through
experimental control, such as gate-voltage tuning 44.

For experimental observations, we propose utilizing scanning magnetometry
techniques with nitrogen-vacancy centers 29, as well as Lorentz transmission
electron microscopy 45 and magnetic transmission soft X-ray microscopy 13, to
directly observe meron pairs ranging in size from 80 to 20 nm at different
twist angles ($\theta$ = 0.5° to 2°). Indirect measurements can involve
detecting anomalous kinks in the magnetization curve, which can serve as an
indication of the presence of merons. Techniques such as the magneto-optical
Kerr effect, commonly employed in the study of 2D vdW magnets 30, 31, can
offer valuable insights for such indirect measurements.

Future research should consider incorporating various magnetic interactions
present in vdW magnets that were overlooked in our current model, such as
exchange anisotropy, interactions beyond nearest neighbors 16, the
Dzyaloshinskii-Moriya interaction, and magnetic dipole-dipole interactions 15.
An important research question to explore is how these additional interactions
impact the stabilization of the meron quartet and their behaviors in vdW
magnets.

We highlight that the discovery and realization of merons in 2D vdW magnets
via twist open a unique avenue for investigating their fascinating properties
with remarkable flexibility, through either external stimuli 35, 36 or the
creation of heterostructures 37. This significant breakthrough not only
deepens our understanding of these fractionalized topological spin textures in
magnets but also holds great promise for future technological advancements.

## References

  *   * 1 Nagaosa, N. & Tokura, Y.  Topological properties and dynamics of magnetic skyrmions.  _Nat. Nanotechnol._ 8, 899–911 (2013).
  * 2 Shinjo, T., Okuno, T., Hassdorf, R., Shigeto, K. & Ono, T.  Magnetic vortex core observation in circular dots of permalloy.  _Science_ 289, 930–932 (2000).
  * 3 Phatak, C., Petford-Long, A. K. & Heinonen, O.  Direct observation of unconventional topological spin structure in coupled magnetic discs.  _Phys. Rev. Lett._ 108, 067205 (2012).
  * 4 Wintz, S. _et al._ Topology and origin of effective spin meron pairs in ferromagnetic multilayer elements.  _Phys. Rev. Lett._ 110, 177201 (2013).
  * 5 Tan, A. _et al._ Topology of spin meron pairs in coupled Ni/Fe/Co/Cu(001) disks.  _Phys. Rev. B_ 94, 014433 (2016).
  * 6 Shigeto, K., Okuno, T., Mibu, K., Shinjo, T. & Ono, T.  Magnetic force microscopy observation of antivortex core with perpendicular magnetization in patterned thin film of permalloy.  _App. Phys. Lett._ 80, 4190–4192 (2002).
  * 7 Fu, X. _et al._ Optical manipulation of magnetic vortices visualized in situ by lorentz electron microscopy.  _Sci. Adv._ 4, eaat3077 (2018).
  * 8 Van Waeyenberge, B. _et al._ Magnetic vortex core reversal by excitation with short bursts of an alternating field.  _Nature_ 444, 461–464 (2006).
  * 9 Ruotolo, A. _et al._ Phase-locking of magnetic vortices mediated by antivortices.  _Nat. Nanotechnol._ 4, 528–532 (2009).
  * 10 Roy, P. E. _et al._ Antivortex domain walls observed in permalloy rings via magnetic force microscopy.  _Phys. Rev. B_ 79, 060407 (2009).
  * 11 Chmiel, F. P. _et al._ Observation of magnetic vortex pairs at room temperature in a planar $\alpha$-Fe2O3/Co heterostructure.  _Nat. Mater._ 17, 581–585 (2018).
  * 12 Yu, X. Z. _et al._ Transformation between meron and skyrmion topological spin textures in a chiral magnet.  _Nature_ 564, 95–98 (2018).
  * 13 Gao, N. _et al._ Creation and annihilation of topological meron pairs in in-plane magnetized films.  _Nat. Commun._ 10, 5603 (2019).
  * 14 Bedoya-Pinto, A. _et al._ Intrinsic 2D-XY ferromagnetism in a van der Waals monolayer.  _Science_ 374, 616–620 (2021).
  * 15 Lu, X., Fei, R., Zhu, L. & Yang, L.  Meron-like topological spin defects in monolayer CrCl3.  _Nat. Commun._ 11, 4724 (2020).
  * 16 Augustin, M., Jenkins, S., Evans, R. F. L., Novoselov, K. S. & Santos, E. J. G.  Properties and dynamics of meron topological spin textures in the two-dimensional magnet CrCl3.  _Nat. Commun._ 12, 185 (2021).
  * 17 Tong, Q., Liu, F., Xiao, J. & Yao, W.  Skyrmions in the moiré of van der Waals 2D magnets.  _Nano Lett._ 18, 7194–7199 (2018).
  * 18 Hejazi, K., Luo, Z.-X. & Balents, L.  Heterobilayer moiré magnets: Moiré skyrmions and commensurate-incommensurate transitions.  _Phys. Rev. B_ 104, L100406 (2021).
  * 19 Akram, M. & Erten, O.  Skyrmions in twisted van der Waals magnets.  _Phys. Rev. B_ 103, L140406 (2021).
  * 20 Xiao, F., Chen, K. & Tong, Q.  Magnetization textures in twisted bilayer CrX3 (X=Br, I).  _Phys. Rev. Research_ 3, 013027 (2021).
  * 21 Akram, M. _et al._ Moiré skyrmions and chiral magnetic phases in twisted CrX3 (X = I, Br, and Cl) bilayers.  _Nano Lett._ 21, 6633–6639 (2021).
  * 22 Ghader, D., Jabakhanji, B. & Stroppa, A.  Whirling interlayer fields as a source of stable topological order in moiré CrI3.  _Commun. Phys._ 5, 192 (2022).
  * 23 Zheng, F.  Magnetic skyrmion lattices in a novel 2D-twisted bilayer magnet.  _Adv. Func. Mater._ 33, 2206923 (2023).
  * 24 Kim, K.-M., Kiem, D. H., Bednik, G., Han, M. J. & Park, M. J.  Ab initio spin hamiltonian and topological noncentrosymmetric magnetism in twisted bilayer CrI3.  _Nano Lett._ 23, 6088–6094 (2023).
  * 25 Fumega, A. O. & Lado, J. L.  Moiré-driven multiferroic order in twisted CrCl3, CrBr3 and CrI3 bilayers.  _2D Materials_ 10, 025026 (2023).
  * 26 Berezinsky, V. L.  Destruction of long range order in one-dimensional and two-dimensional systems having a continuous symmetry group. I. Classical systems.  _Sov. Phys. JETP_ 32, 493–500 (1971).
  * 27 Kosterlitz, J. M. & Thouless, D. J.  Ordering, metastability and phase transitions in two-dimensional systems.  _J. Phys. C_ 6, 1181–1203 (1973).
  * 28 Burch, K. S., Mandrus, D. & Park, J.-G.  Magnetism in two-dimensional van der Waals materials.  _Nature_ 563, 47–52 (2018).
  * 29 Song, T. _et al._ Direct visualization of magnetic domains and moiré magnetism in twisted 2D magnets.  _Science_ 374, 1140–1144 (2021).
  * 30 Xu, Y. _et al._ Coexisting ferromagnetic–antiferromagnetic state in twisted bilayer CrI3.  _Nat. Nanotechnol._ 17, 143–147 (2022).
  * 31 Xie, H. _et al._ Twist engineering of the two-dimensional magnetism in double bilayer chromium triiodide homostructures.  _Nat. Phys._ 18, 30–36 (2022).
  * 32 Xie, H. _et al._ Evidence of non-collinear spin texture in magnetic moiré superlattices.  _Nat. Phys._ (2023).
  * 33 Hejazi, K., Luo, Z.-X. & Balents, L.  Noncollinear phases in moiré magnets.  _PNAS_ 117, 10721–10726 (2020).
  * 34 Kim, K.-M. & Park, M. J.  Controllable magnetic domains in twisted trilayer magnets.  Preprint at arXiv https://doi.org/10.48550/arXiv.2306.14423 (2023).
  * 35 Wang, Z. _et al._ Electric-field control of magnetism in a few-layered van der Waals ferromagnetic semiconductor.  _Nat. Nanotechnol._ 13, 554–559 (2018).
  * 36 Jiang, S., Li, L., Wang, Z., Mak, K. F. & Shan, J.  Controlling magnetism in 2D CrI3 by electrostatic doping.  _Nat. Nanotechnol._ 13, 549–553 (2018).
  * 37 Gibertini, M., Koperski, M., Morpurgo, A. F. & Novoselov, K. S.  Magnetic 2D materials and heterostructures.  _Nat. Nanotechnol._ 14, 408–419 (2019).
  * 38 Sivadas, N., Okamoto, S., Xu, X., Fennie, C. J. & Xiao, D.  Stacking-dependent magnetism in bilayer CrI3.  _Nano Lett._ 18, 7658–7664 (2018).
  * 39 Jang, S. W., Jeong, M. Y., Yoon, H., Ryee, S. & Han, M. J.  Microscopic understanding of magnetic interactions in bilayer CrI3.  _Phys. Rev. Mater._ 3, 031001 (2019).
  * 40 Song, T. _et al._ Switching 2D magnetic states via pressure tuning of layer stacking.  _Nat. Mater._ 18, 1298–1302 (2019).
  * 41 Li, T. _et al._ Pressure-controlled interlayer magnetism in atomically thin CrI3.  _Nat. Mater._ 18, 1303–1308 (2019).
  * 42 Chen, W. _et al._ Direct observation of van der Waals stacking–dependent interlayer magnetism.  _Science_ 366, 983–987 (2019).
  * 43 Hubert, A. & Schäfer, R.  _Magnetic Domains_ (Springer Berlin Heidelberg, 1998).
  * 44 Tang, M. _et al._ Continuous manipulation of magnetic anisotropy in a van der Waals ferromagnet via electrical gating.  _Nat. Electron._ 6, 28–36 (2023).
  * 45 Ding, B. _et al._ Observation of magnetic skyrmion bubbles in a van der Waals ferromagnet Fe3GeTe2.  _Nano Lett._ 20, 868–873 (2020).

##  4 Methods

###  4.1 Moiré superlattice

We utilized commensurate moiré superlattice structures generated by rotating
one layer of an aligned honeycomb lattice bilayer relative to the other layer.
The rotation is centered at one of the hexagonal centers and is defined by the
angle
$\theta=\arccos{\left(\frac{m^{2}+n^{2}+4mn}{2m^{2}+2n^{2}+2mn}\right)}$. Our
specific focus was on the case of $n=m+1$, where the superlattice periodicity
aligns with the periodicity of moiré patterns.

###  4.2 Interlayer Heisenberg exchange interactions

Figure 6: Illustration of the procedure to determine interlayer Heisenberg
exchange interactions. a-f: Each panel depicts the procedure for determining
interlayer Heisenberg exchange interactions $J_{ij}^{\perp}$ in different
local stacking patterns. The left images show the local stacking patterns,
with red and blue spheres representing the positions of local spins in the top
and bottom layers, respectively. Yellow (blue) arrows indicate the parallel
coordinate displacements from A (B) sublattice spins to their neighboring
spins. The middle (right) color plot labeled “Same sublattice” (“Different
sublattice”) depicts $J_{ij}^{\perp}$ as a function of the parallel coordinate
displacements between two spins on the top and bottom layers
$(x_{ij},y_{ij})=(x_{i}-x_{j},y_{i}-y_{j})$, where red and blue colors
indicate AFM and FM characters, respectively. The same displacement vectors
shown in the left images are drawn onto these maps, illustrating the
interpolation process for determining $J_{ij}^{\perp}$ corresponding to each
bond. The coupling maps were adopted from previous ab initio calculations on
bilayer CrI3 24.

We determined $J_{ij}^{\perp}$ in our spin model by utilizing the coupling
maps for interlayer Heisenberg exchange interactions derived from previous ab
initio calculations on bilayer CrI3 24, as depicted in the color plots of Fig.
6. Our method 23, 24 involves computing the parallel coordinate displacements
$(x_{ij},y_{ij})$ for all bonds between spins from the top and bottom layers,
respectively. These displacements are then interpolated onto the coupling
maps, which are functions of $(x_{ij},y_{ij})$, to obtain the values of
$J_{ij}^{\perp}$ for each bond. Fig. 6a-f illustrate this procedure for
representative stacking patterns within the moiré superlattice. There are two
distinct types of bonds: (i) bonds between two spins on the same sublattice
and (ii) bonds between two spins on different sublattices. The coupling maps
labeled “Same sublattice” and “Different sublattice” are utilized for the
former and latter cases, respectively. The interactions for bond distances
$\sqrt{x_{ij}^{2}+y_{ij}^{2}}$ larger than $\sqrt{3}a\approx 5.2$Å are
neglected, as their contribution is negligible.

###  4.3 Iterative optimization method

Figure 7: Relaxation of spin configurations using the iterative optimization
method. Panel a represents the MD state, as depicted in Fig. 4d-f. Panel b
corresponds to the MQ state, as illustrated in Fig. 4a-c. The blue curves
represent the total magnetic energy, and the red curves represent the
convergence factor in the relaxation defined as
$\frac{1}{N}\sum_{i}|\mathbf{h}_{i}^{\perp}|$.

We obtained a stable spin configuration using an iterative optimization
method. In this procedure, we update a trial configuration by subtracting the
normalized spin vector $\mathbf{n}_{i}=\mathbf{S}_{i}/|\mathbf{S}_{i}|$ with
the gradient vector of the spin Hamiltonian $H$, denoted as
$\mathbf{h}_{i}\equiv\frac{\partial
H[\\{\mathbf{n}_{i}\\}]}{\partial\mathbf{n}_{i}}$. Here, $\mathbf{S}_{i}$
represents the classical spin vector at the $i$-th site. The update rule is
given by:

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

where $(s)$ denotes the iteration step, and $c$ controls the update rate. The
iterative process continues until $\mathbf{h}_{i}^{(s)}$ becomes parallel to
$\mathbf{n}_{i}^{(s)}$, implying that further updates will not change the spin
configuration. At this stage, the magnetic energy is minimized. To ensure
numerical convergence, we used the condition
$\frac{1}{N}\sum_{i}|\mathbf{h}_{i}^{\perp}|<10^{-10}$, where
$\mathbf{h}_{i}^{\perp}=\mathbf{h}_{i}-\mathbf{n}_{i}(\mathbf{n}_{i}\cdot\mathbf{h}_{i})$.
This convergence condition is identical to the stability condition in the
Landau-Lifshitz-Gilbert equation describing the magnetic dynamics. We
incorporated periodic boundary conditions during the iterative updates to
maintain consistency across the boundaries of the moiré unit cell.

###  4.4 Local magnetic energy maps

We computed the local energy maps in Fig. 4 and Supplementary Videos (5.1,
5.2) by utilizing the following equations:

$\displaystyle E_{\textrm{tot}}(\mathbf{r}_{i})$
$\displaystyle=E_{\textrm{t}}(\mathbf{r}_{i})+E_{\textrm{b}}(\mathbf{r}_{i})+E_{\textrm{tb}}(\mathbf{r}_{i}),$
(5a) $\displaystyle E_{\textrm{t}}(\mathbf{r}_{i})$
$\displaystyle=-\frac{JS^{2}}{2}\sum_{j\in\textrm{N.N}(i)}\mathbf{n}_{i}^{t}\cdot\mathbf{n}_{j}^{t}+AS^{2}\left(\mathbf{n}_{i}^{t}\cdot\hat{z}\right)^{2},$
(5b) $\displaystyle E_{\textrm{b}}(\mathbf{r}_{i})$
$\displaystyle=-\frac{JS^{2}}{2}\sum_{j\in\textrm{N.N}(i)}\mathbf{n}_{i}^{b}\cdot\mathbf{n}_{j}^{b}+AS^{2}\left(\mathbf{n}_{i}^{b}\cdot\hat{z}\right)^{2},$
(5c) $\displaystyle E_{\textrm{tb}}(\mathbf{r}_{i})$
$\displaystyle=S^{2}\sum_{j}J_{ij}^{\perp}\mathbf{n}_{i}^{t}\cdot\mathbf{n}_{j}^{b}.$
(5d)

Here, $\mathbf{r}_{i}$ represents the position of $\mathbf{n}_{i}$, and
$\text{N.N}(i)$ denotes the nearest neighbors of the $i$-th site. We have set
$S=1$ for simplicity.

###  4.5 Determination of FM-MD phase boundary

Figure 8: Ferromagnetic order parameter ($M$), defined by Eq. (6), plotted as
a function of twist angle ($\theta$) for various intralayer exchange values
($J$) (in meV unit).

We utilized random initial configurations represented as
$\\{\mathbf{n}_{i}^{(0)}\\}=\left\\{(u_{i},v_{i},w_{i})\right\\}$, where
$u_{i}$, $v_{i}$, and $w_{i}$ are random numbers within the range of $[-1,1]$.
Normalization was performed to ensure unit length for all
$\mathbf{n}_{i}^{(0)}$. These initial configurations were relaxed through the
iterative optimization method to obtain the ground state configuration
$\\{\mathbf{n}_{i}\\}$ for each $\theta$ and $J$. From these ground states, we
computed the ferromagnetic order parameter defined as:

$\displaystyle M=\frac{1}{N}\left|\sum_{i}\mathbf{n}_{i}\right|.$ (6)

The results are shown in Fig. 8. From these results, we identified a threshold
value of $\theta$ for each $J$ such that $M<1$. These values were utilized to
determine the phase boundary between the FM and MD phases, as depicted in Fig.
2c of the main text.

###  4.6 Critical twist angle formula for FM-MD transition

We derived the critical twist angle formula using an effective energy
functional in the continuum limit, corresponding to the spin Hamiltonian in
Eq. (2.1) of the main text:

$\displaystyle E[\mathbf{n}_{l}]$ $\displaystyle=\frac{1}{2}J\sum_{l=t,b}\int
d^{2}\mathbf{r}\left|\nabla\mathbf{n}_{l}(\mathbf{r})\right|^{2}+\frac{A}{a^{2}}\sum_{l=t,b}\int
d^{2}\mathbf{r}\left(n_{l,z}(\mathbf{r})\right)^{2}$
$\displaystyle+\frac{1}{a^{2}}\int
d^{2}\mathbf{r}~{}J_{\perp}\left(\mathbf{r}\right)\mathbf{n}_{t}(\mathbf{r})\cdot\mathbf{n}_{b}(\mathbf{r}),$
(7)

where $\mathbf{n}_{l}(\mathbf{r})$ represents the vector field depicting the
orientation of spins at position $\mathbf{r}$ in the top ($l=t$) and bottom
($l=b$) layers, respectively. To compute the magnetic energy for the MD phase,
we employed several approximations. Firstly, we considered the shape of an AFM
patch as a disk with a size of $L$. Secondly, we employed the following
function form for the interlayer coupling $J_{\perp}(\mathbf{r})$:

$\displaystyle
J_{\perp}(\mathbf{r})=J_{\perp}^{\textrm{max}}(1-|\mathbf{r}|/L),$ (8)

where $J_{\perp}^{\textrm{max}}$ is the maximum value of
$J_{\perp}(\mathbf{r})$ at the center of the patch. Lastly, we assumed
$\mathbf{n}_{t,b}(\mathbf{r})=(0,\sin{\phi_{t,b}(\mathbf{r})},\cos{\phi_{t,b}(\mathbf{r})})$
with the following function forms for the phase angles
$\phi_{t,b}(\mathbf{r})$:

$\displaystyle\phi_{t}(\mathbf{r})$
$\displaystyle=+\phi_{0}(1-|\mathbf{r}|/L),$
$\displaystyle\phi_{b}(\mathbf{r})$
$\displaystyle=-\phi_{0}(1-|\mathbf{r}|/L),$ (9)

where $\phi_{t,b}(\mathbf{r})$ reach their maximum values at the center of the
disk. Inserting Eqs. (8) and (4.6) into Eq. (4.6), we obtained an energy
function as a function of $\phi_{0}$ as

$\displaystyle E(\phi_{0})$ $\displaystyle=\pi
J\phi_{0}^{2}+\frac{\bar{J}_{\perp}\pi
L^{2}}{a^{2}}\bigg{(}\frac{3\cos\phi_{0}(\sin\phi_{0}-\phi_{0}\cos\phi_{0})}{\phi_{0}^{3}}\bigg{)}$
$\displaystyle\approx\frac{\bar{J}_{\perp}\pi
L^{2}}{a^{2}}+\phi_{0}^{2}\bigg{(}\pi J-\frac{3\bar{J}_{\perp}\pi
L^{2}}{5a^{2}}\bigg{)}+\mathcal{O}(\phi_{0}^{4}),$ (10)

where $\bar{J}_{\perp}$ is the average value of defined as
$\bar{J}_{\perp}=\frac{2\pi}{\pi
L^{2}}\int^{L}_{0}drrJ_{\perp}(\mathbf{r})=J_{\perp}^{\textrm{max}}/3$. When
the sign of the quadratic term changes, the phase transition occurs. Equating
the coefficient of the quadratic term to zero, we find a critical AFM patch
size $L_{c}$:

$\displaystyle L_{c}=a\sqrt{\frac{5J}{3\bar{J}_{\perp}}}.$ (11)

From the relation $3\pi L^{2}=f\frac{\sqrt{3}}{2}L_{\textrm{M}}^{2}$, where
$f$ is the fraction of spins in the AFM domains and
$L_{\textrm{M}}=\frac{a}{2\sin(\theta/2)}\approx\frac{a}{\theta}$ is the
lattice constant of the moiré superlattice ($3$ comes from the fact that there
are three AFM patches per moiré unit cell), we find the relation
$L=\sqrt{\frac{f\sqrt{3}}{6\pi}}L_{\textrm{M}}\approx\sqrt{\frac{f\sqrt{3}}{6\pi}}\frac{a}{\theta}$.
Using this relation with Eq. (11), we finally obtain the critical twist angle
$\theta_{c}$ as

$\displaystyle\theta_{c}=\sqrt{\frac{\sqrt{3}}{10\pi}\frac{f\bar{J}_{\perp}}{J}}.$
(12)

We utilized this formula to fit the numerically calculated FM-MD phase
boundary, as presented in Fig. 2c of the main text. The fitting parameter was
obtained as $f\bar{J}_{\perp}\approx 0.56$meV.

###  4.7 Determination of MD-MQ phase boundary

Figure 9: Left: Ferromagnetic order parameter ($M$) and Right: total meron
number ($N_{\textrm{meron}}$), as defined in Eqs. (6) and (13), computed from
different stable magnetic states. Each state was obtained from a distinct
random initial spin configuration. The parameter values of $J=2$meV and
$\theta=1.61$° are utilized.

We utilized the same relaxed configurations that were employed to determine
the FM-MD phase boundary, as depicted in Sec. 4.5. For each spin
configuration, we computed the total number of merons from both the top and
bottom layers by using the following formula:

$\displaystyle N_{\textrm{meron}}=2\times\sum_{l=t,b}\frac{1}{4\pi}\int dx\int
dy\left|(\partial_{x}\mathbf{n}_{l}\times\partial_{y}\mathbf{n}_{l})\cdot\mathbf{n}_{l}\right|,$
(13)

where the integral is performed throughout the moiré unit cell. The absolute
value was imposed in order to consider both the number of merons and the
number of antimerons, preventing their skyrmion numbers from canceling out
with each other. The factor of $2$ comes from the fact that the contribution
of a meron or an antimeron is $1/2$. The results for specific parameter values
of $J=2$ meV and $\theta=1.61$° are displayed in Fig. 9. We confirmed that the
value of $N_{\textrm{meron}}\approx 4$ corresponds to the actual formation of
the MQ state. By establishing this correspondence, we computed the total meron
number in each relaxed configuration for different values of $J$ and $\theta$.
For each $J$, we identified a threshold value of $\theta$ for the appearance
of $N_{\textrm{meron}}\approx 4$. These results were utilized to determine the
phase boundary between the MD and MQ phases, as depicted in Fig. 5a of the
main text.

###  4.8 Determination of critical field strengths for MQ-MD transition

Figure 10: Effect of the Zeeman term, as defined in Eq. (3) of the main text,
on the modification of the MQ state, as depicted in Fig. 4a-c. The left and
right axes correspond to the total meron number ($N_{\textrm{meron}}$) and
magnetic energy, as defined in Eqs. (13) and (5a), respectively. The parameter
$B$ represents an external magnetic field in the out-of-plane direction.
$B_{c1}$ and $B_{c2}$ denote the critical field strengths for the transition
of MQ state to the MD state. The parameter values of $J=2$meV and
$\theta=1.61$° are utilized.

In the first step, we utilized random initial configurations represented as
$\\{\mathbf{n}_{i}^{(0)}\\}=\\{(0,0,1)+0.1*(u_{i},v_{i},w_{i})\\}$, where
$u_{i}$, $v_{i}$, and $w_{i}$ are random numbers within the range of $[-1,1]$.
Normalization was performed to ensure unit length for all
$\mathbf{n}_{i}^{(0)}$. The utilization of these polarized configurations
ensured consistency in our analysis by imposing that all MQ states belong to
the same class with the same polarity $p=-1/2$ among the meron cores. These
initial configurations were relaxed through the iterative optimization method.
Through this procedure, we prepared more than 20 different samples of the MQ
state in the zero-field condition.

In the subsequent step, we proceeded to compute modified spin configurations
from these samples by incorporating the Zeeman term described in Eq. (3) of
the main text. We adjusted the field strength from $B=-0.86T$ to $B=0.86T$ in
this term, as depicted in Fig. 5b of the main text. We used $g=2$ and
$\mu_{B}=0.058$meV, and for simplicity, we assumed $S=1$. Each sample was
relaxed under the influence of the Zeeman term, using an iterative
optimization method. After the relaxation, we computed the magnetic energy,
the total number of merons, and the net magnetization in the out-of-plane
direction, using Eqs. (5a) and (13), and $M_{z}=\frac{1}{N}\sum_{i}n_{i,z}$,
respectively. We repeated this computation for different field strengths in
order to examine the effects of varying magnetic fields. Fig. 10 illustrates
the results for the magnetic energy and total number of merons. The physical
quantities consistently exhibit discontinuous jumps at critical field
strengths denoted by $B_{c1}$ and $B_{c2}$. We confirmed the discontinuity
corresponds to the actual transition from the MQ to the MD states (we checked
that the small jump before $B_{c2}$ corresponds to the reversal of the
polarity of the merons). Utilizing this correspondence, we computed the
critical field strength for each sample and each parameter of $J$ and
$\theta$. These were utilized to determine the critical field strengths for
the MQ-MD transition, as depicted in Fig. 5c of the main text.

###  4.9 MQ states with different skyrmion numbers

Figure 11: Illustrative examples of MQ states showcasing different skyrmion
numbers $(Q_{t},Q_{b})$ in the top and bottom layers, as defined in Eq. (14).
a-b: $(Q_{t},Q_{b})=(1,1)$, c-d: $(Q_{t},Q_{b})=(1,0)$, and e-f:
$(Q_{t},Q_{b})=(1,-1)$. Arrows indicate in-plane components, while the color
scale in the markers represents out-of-plane components. Marker sizes are
adjusted for better visibility. Shaded areas indicate AFM patches. The
parameters $J=2$ meV and $\theta=1.61$° are utilized.

Fig. 11 presents illustrative examples of MQ states, which exhibit distinct
values of skyrmion numbers $(Q_{t},Q_{b})$ in the top and bottom layers,
respectively, defined as,

$\displaystyle Q_{t}$ $\displaystyle=\frac{1}{4\pi}\int dx\int
dy(\partial_{x}\mathbf{n}_{t}\times\partial_{y}\mathbf{n}_{t})\cdot\mathbf{n}_{t},$
(14a) $\displaystyle Q_{b}$ $\displaystyle=\frac{1}{4\pi}\int dx\int
dy(\partial_{x}\mathbf{n}_{b}\times\partial_{y}\mathbf{n}_{b})\cdot\mathbf{n}_{b}.$
(14b)

The first row (a-b) illustrates a pair of two merons in the top layer
($\textrm{M}_{\textrm{t1}}$–$\textrm{M}_{\textrm{t1}}$) and a pair of two
antimerons
($\overline{\textrm{M}}_{\textrm{b1}}$–$\overline{\textrm{M}}_{\textrm{b2}}$)
in the bottom layer, corresponding to $(Q_{t},Q_{b})=(1,-1)$. The second row
(c-d) shows a pair of two merons in the top layer
($\textrm{M}_{\textrm{t1}}$–$\textrm{M}_{\textrm{t2}}$) and a meron-antimeron
pair ($\textrm{M}_{\textrm{b}}$–$\overline{\textrm{M}}_{\textrm{b}}$) in the
bottom layer, corresponding to $(Q_{t},Q_{b})=(1,-1)$. Lastly, the third row
(e-f) exhibits two pairs of two antimerons in the top layer
($\overline{\textrm{M}}_{\textrm{t1}}$–$\overline{\textrm{M}}_{\textrm{t2}}$)
and
($\overline{\textrm{M}}_{\textrm{b1}}$–$\overline{\textrm{M}}_{\textrm{b2}}$)
in the bottom layer, corresponding to $(Q_{t},Q_{b})=(1,-1)$. The MQ state in
Fig. 3c-d of the main text features $(Q_{t},Q_{b})=(0,0)$.

##  5 Supplementary information

###  5.1 Supplementary Video 1

Relaxation process for the MQ state, as presented in Fig. 4a-c and g-i. First
row: Spin configurations in the top and bottom layers, and relative
orientation. Text box in the top right corner displaying the iteration step
and the magnetic energy per spin at each step. Second row: Intralayer exchange
energy plus single-ion anisotropy energy in the top and bottom layers, and
interlayer exchange energy.

###  5.2 Supplementary Video 2

Relaxation process for the MD state, as presented in Fig. 4d-f and j-l. First
row: Spin configurations in the top and bottom layers, and relative
orientation. Text box in the top right corner displaying the iteration step
and the magnetic energy per spin at each step. Second row: Intralayer exchange
energy plus single-ion anisotropy energy in the top and bottom layers, and
interlayer exchange energy.

Data availability All data generated or analysed during this study are
included in this Article. Source data are provided with this paper.

Code availability The code written for use in this study is available from the
corresponding author upon reasonable request.

Acknowledgements

K.K. acknowledges support from the Institute for Basic Science in the Republic
of Korea under the project IBS-R024-D1. G.G. was supported by the National
Research Foundation of Korea (NRF-2022R1C1C2006578). M.J.P. was supported by
the National Research Foundation of Korea (NRF) grant funded by the Korea
government (MSIT) (RS-2023- 00252085,RS-2023-00218998). Lastly, S.K.K. was
supported by Brain Pool Plus Program through the National Research Foundation
of Korea funded by the Ministry of Science and ICT (2020H1D3A2A03099291) and
National Research Foundation of Korea funded by the Korea Government via the
SRC Center for Quantum Coherence in Condensed Matter (RS-2023-00207732).


# Emergence of autopoietic vesicles able to grow, repair and reproduce in a
minimalist particle system

Thomas Cabaret
Nice, France
Email<EMAIL_ADDRESS>
###### Abstract

This paper describes a 2D particle system in which autopoietic vesicles that
are able to grow, repair, and reproduce can emerge from an initial random
distribution of individual components. The reproduction process also exhibits
inheritance of the vesicle content that could in principle change from one
generation to the next and influence fitness, making those autopoietic
vesicles eligible to be acted upon by natural selection. Video abstract: link.

## Keywords

Origin of life, Simulation, Particle system, Autopoiesis, Autocatalytic set,
Abiogenesis, A-Life.

##  1 Introduction

Using artificial environments to witness some form of "origin of life" events
has been the subject or prior works in the cellular automata world[1][2], or
with particle systems like PPS[3] or Squirm3[4]. [1] is a variation of the
Langton’s loops showing spontaneous emergence of replicators from an initial
random distribution, and a limited form of evolution. It comes with highly
handcrafted rules and bypasses important real world constraints like
conservation of matter (conservation of the number of particles). Evoloop[2]
is another variation of Langton’s loop that stresses a richer evolvability, at
the expense of losing spontaneous formation capability. On the particle system
side, PPS[3] highlights the relationship between rulesets and emergence of
replicators. While it aligns more with real-world constraints, the
spontaneously formed structures currently lack a clear mechanism for evolution
and extending the system towards this goal would be challenging. Conversely,
Squirm3[4], another particle system, exhibits highly ordered structures able
to reproduce with information inheritance from one generation to the next, but
those structures cannot form spontaneously. Here we describe a 2D particle
system in which the rules respect conservation of matter and lead to the
spontaneous formation of vesicular structures able to grow, repair, and
reproduce. Our particle system tries to find the right balance between fine-
tuning and emergence and focuses on the stability and higher order properties
of those emergent replicators, like their ability to convey some form of
information from one generation to the next and so to evolve.

##  2 Definitions

Prior to describing our particle system and corresponding findings, it is
essential to introduce two fundamental concepts that underpin its
construction, autocatalytic sets and autopoiesis, which play important roles
in various other origin of life research.

We will refer to any set of chemical reactants as an autocatalytic set if, for
any reactants (catalysts included) not provided by the environment, there is a
reaction within the set that produces those reactants and ensures their
concentration increases over time.

Figure 1: Living system as an autopoietic structure, adapted from [5] Figure
2: a), b), c) Simulation 1: starting from a random distribution of individual
components and leading to non empty vesicle formation. d), e), f) Simulation
2: starting from a fully formed autopoietic vesicle exposed to a source of
precursor particles and leading to the vesicle reproduction.

A structure will be said autopoietic when it is able to produce its components
and repair itself from within. See figure 1 from (Luisi 2003)[5] where a
living system can be understood as the alliance of structure and inner
workings, such that the structure enables the inner workings, and the inner
workings, in turn, sustain the structure. In addition of using autocatalysis,
autopoietic structures are also able to maintain their organization and
integrity. It can also be noted that autopoiesis has an implicit relationship
with reproduction, as it can allow growth, division, and repair of daughter
parts thus achieving a full reproduction cycle.

##  3 Simulation and findings

In order to gain insight into the particle system we have developed and the
significant dynamics exhibited in various simulations, we suggest to start by
examining one such simulation and provide a cursory explanation of the
different rules governing our particle system as we progress (This simulation
is available as a YouTube video see annex 2). A more exhaustive and technical
presentation of these rules will be reserved for the next chapter.

We will start a simulation from a random distribution of particles and
describe the various stages of chemical evolution, in addition to highlighting
the characteristics of autonomously formed structures able to grow, repair,
and reproduce.

| Precursor | Catalyst | Intermediary | Surfactant
---|---|---|---|---
Precursor | none | only chemical | only chemical | none
Catalyst | only chemical | none | none | close range repulsion
Intermediary | only chemical | none | none | close range repulsion
Surfactant | none | close range repulsion | close range repulsion | complex force profile, see text for details
Table 1: Summary of interactions between particle types.

In figure 2a we see such an initial random distribution of the 4 types of
particles used in our model: green particles will be referred to as
Precursors, red particles will be referred to as Catalysts, purple particles
as Intermediary, and yellow particles as Surfactant or polar particles. After
a short while, like illustrated by figure 2b, we see small aggregates of
yellow particles that exhibit an intrinsic curvature. On figure 2c we see that
when enough surfactant particles aggregate full vesicles are eventually
formed, some of which can entrap either at least a catalyst or an intermediary
particle. A second simulation, starting from a fully formed vesicle, shows
that those vesicles containing catalysts are indeed autopoietic structures
able to grow, repair, and reproduce. When such a vesicle is placed close to a
source of precursor particles, those can, by design, cross the surfactant
membrane, and a proto-metabolic activity can start inside, allowed by a set of
chemical rules changing types of particles when specific ones come close to
each other. The net effect of the inner chemical reactions is to consume
precursors and increase the inner population of catalyst, intermediary and
surfactants particles that eventually aggregate to the membrane from within,
making the whole structure grow. See figure 2d. Large vesicles are unstable,
start to squirm and eventually break apart into 2 halves that initially stay
close to each other, as illustrated in figure 2e. This proximity usually
prevents the leaking of a significant part of the content. An autopoietic
activity persists in both halves, producing surfactants that leads quickly to
both vesicle repair, achieving a full reproduction cycle. See figure 2f.

##  4 Model and genesis

Our particle system involves two chemical rules, see figure 3a: A Precursor
can be transmuted into an Intermediary particle in presence of a Catalyst, and
if a Precursor and an Intermediary particle interacts they are both transmuted
into Catalyst and Surfactant. This set of chemical rules allow the existence
of a minimalist autocatalytic set, see figure 3b. In the simulation both of
those chemical rules are triggered if eligible particles are within an
interaction radius, small enough to prevent most of the chemical interactions
between particles located on opposite sides of a surfactant membrane.

Figure 3: a) The two chemical rules of our particle system. b) The minimalist
autocatalytic set they define. Bold arrows for chemical transformations,
dashed arrows for catalytic activity and gray arrows for external intake. One
color dot by particle type: precursor(green), catalyst(red),
intermediary(purple) and surfactant(yellow).

Our particle system also involves forces between particles that can modify
their trajectories. The force profile between particles is summarized in table
1. The only non trivial force profile concerns surfactant molecules, it’s
carefully handcrafted to lead to a satisfactory vesicular behavior. The key
idea to understand this force profile is to think of it as derived from a
potential with minimum being reached when surfactant molecules are close to
each other, side by side, tail on the same side, almost parallel to each other
with a slight angle. See figure 4a.

Figure 4: Surfactant molecules force profile details a) equilibrium state, b)
repulsive domain.

If a couple of interacting yellow particles are far from this equilibrium
state, each will experience a force towards it. The slight deviation from
perfect parallelism in this stable equilibrium state (figure 4a) is a way to
impose intrinsic local curvature to surfactant molecule aggregates. This
deviation being constant, too large vesicles are by nature unstable as all of
their surfactant molecules cannot be at equilibrium with their neighbors
simultaneously. It was designed that way on purpose with the hope it would
lead to some kind of division and the first behavior observed was
satisfactory, as seen previously, large vesicles break apart and then repair.
The force profile of surfactant molecules also contain a repulsive domain when
they are close, head to head, and anti-parallel, see figure 4b. The key idea
of this specificity is to avoid fusion of vesicles freshly formed as at the
interaction location between two vesicles, surfactant molecules will be in
this repulsive domain, making fully formed vesicles repulsive to one another
which strongly favors the behavior towards division rather than fusion.

Figure 5: Visualization of force interaction radius illustrated by the black
circles around surfactant particles.

All the force interaction radius are small, as illustrated by the circles in
figure 5 showing interaction radius of surfactant particles which barely allow
each of them to see the neighbor of their neighbor within a vesicle. It is
particularly noteworthy that the interaction range is substantially smaller
than the size of the structures of interest, namely the vesicles. This
significant difference in scale between local interactions and the overall
structure justifies to talk about emergence when referring to the vesicles and
their dynamic.

We can then conclude the description of our particle system with additional
remarks

  * •

Energy is not conserved: the model contains dissipative terms that are
unavoidable to observe structure formation without simulating solvent
particles. The model also includes the exact opposite, Brownian random noise
to favor macroscopic motion of particles and thus particle encounter.

  * •

Momentum is conserved (up to numerical approximations).

  * •

Angular momentum not conserved: the force profile of surfactant molecules can
in some cases lead them to spiral one toward the other to end up side by side.

For more details about the model, please refer to the c++ code in annex.

##  5 Discussion

Despite not having gathered specific metric on the following points, various
manually started simulations show a high stability of vesicle formation from
an initial random distribution: with an initial state involving several
hundreds of particles it is common to see half a dozen fully formed and non
empty vesicles thus making this period of chemical evolution highly
reproducible.

The reproduction seems also quite resilient, although it can sometimes fail,
by example through leakage of the whole content, success seems more frequent
than failure, and most importantly it almost never happens that a reproduction
failure leads to the loss of both daughter vesicles, most of the times at
least one remains from which a reproduction can be retried. In average the
reproduction success rate seems compatible with an exponential growth of the
vesicle population, it is quite common to witness several successful divisions
happening sequentially.

A limitation of our particle system is the need to split the simulation into
two distinct ones: the first for vesicles formation, the second for vesicles
reproduction. This could in principle be fixed by allowing larger simulations
with local macroscopic variation of the environment, like it is the case in
real life by example with a pool at some location and the ocean at another. We
could also introduce conditions that vary periodically with time to mimic real
life cycles. In the current state of the model, if a source of precursor
particles would be added into the initial random starting distribution, those
precursors would mostly react outside emergent vesicles before reaching the
interior of one of them and allowing a more interesting activity.

During various experiments it appeared quite clear that the key for a lot of
relevant behaviors is the existence of the membrane itself. Having particles
able to aggregate into vesicles with a flexible enough dynamic is essential to
define macroscopic systems, to unite unrelated entities into a single one.
Although membrane formation can emerge from numerous sets of rules, stepping
among one of them by pure luck and still have a computationally efficient
model seems unlikely to us. We suggest to study this kind of particle system,
in which vesicular behavior is manually crafted to study higher order
properties, as the dedicated class of vesicular particle systems, which seems
to be justified by the unique macroscopic dynamic they allow.

Finally our simulation could be chemically improved by introducing other types
of particles allowing in principle the evolution of our vesicles as it is
discussed in the next chapter.

##  6 Evolvability of Autopoietic Vesicles

A key question that can arise having structure able to grow and reproduce
concerns Darwinism: can such structures be acted upon by natural selection ?
To show that it is the case we have to show that

  * •

Some information is inherited.

  * •

This information can undergo variations from one generation to another.

  * •

Such variations can influence fitness

The inheritance mechanism of our vesicles is quite straightforward: their
content is conveyed from one generation to the next: the set of types of
particles is conserved with a high probability, and to a lesser extent the
relative concentration between types of particles. If variations of the
vesicle content that are persistent from one generation to the next can occur,
it is clear that this could lead to a variation of the inner chemical activity
and thus impact fitness through a huge number of mechanisms as is observed in
real-world living cells. The current version of our simulation only uses 4
types of particles, but it is in theory possible, by just making the chemistry
rules of the simulation richer, to allow at least one of two following
mechanisms for variation.

###  6.1 Autocatalytic extension

Variation of the content can be achieved by serendipitous leakage within a
vesicle. The new content can become part to the information that is
transferred from one generation to the next if as a whole it forms a new
autocatalytic set. We suggest to refer to this as "autocatalytic extension".
see figure 6.

Figure 6: Example of change by autocatalytic extension: addition of new x, y
types of particles within a vesicle can lead to a new autocatalytic set able
to sustain itself from one generation to the next. Bold arrows for chemical
transformations, dashed arrows for catalytic activity and gray arrows for
external intake. One color dot by particle type: precursor(green),
catalyst(red), intermediary(purple) and surfactant(yellow).

NB: during some simulations we noted an interesting mechanism by which new
particles can enter a vesicle, they can be trapped between the vesicle and a
small chunk of membrane that aggregate to it as illustrated in figure 7.

###  6.2 Programmable reactants

Another way to achieve variation is the introduction of programmable
molecules. The particles of our simulation are more in analogy to real world
molecules than they are to just atoms, and so, some of those types of
particles could be polymers made of several types of elements and able to
catalyze their own formation by a template-copy mechanism. See figure 8. The
introduction of those types of structures within the chemistry would open
entities to an evolution of evolvability as, prior to that, only the catalytic
extension mechanism is a way to mutate and rely on occurrence of rare events
depending on the environment. The introduction of programmable polymer within
the proto-metabolism would allow copy errors from within, thus increasing the
opportunity for variation without relying on external factors.

A suitable framework for implementing the aforementioned enhancements to the
model could be Christian Heinemann’s Alien-project.

Figure 7: Mechanism for vesicle entry discovered during simulation
(illustration): outer particles can end up within vesicle by being trapped
between membrane chunks during their fusion with vesicles. Figure 8: Example
of change by error in template copy. In this figure, the catalyst in red is a
macroscopic abstraction of a polymer made of green and purple particles. It
catalyses it’s own formation by canonical template copy. Bold arrows for
chemical transformations, dashed arrows for catalytic activity and gray arrows
for external intake. One color dot by particle type: precursor(green),
catalyst(red), intermediary(purple) and surfactant(yellow).

##  7 Conclusion

In this paper, we have presented a 2D particle system that exhibits emergent
formation of autopoietic vesicles capable of growth, repair, and reproduction.
We have shown that it is possible for structures with these properties to
emerge from a simple set of rules governing particle interactions. We
highlighted that these structures exhibit inheritance of their content and
discussed how the system could be extended to allow relevant changes from one
generation to the next, thus allowing those structure to be acted upon by
natural selection. While our simulation is highly simplified compared to real-
world chemical systems, it serves as a proof-of-concept that such behaviors
are possible by giving an new A-life example of a minimalist origin of life
event. Our model provides a foundation for further exploration of these
concepts and could be extended and refined to better mimic complex chemistry
and to explore other potential mechanisms for the emergence and evolution of
autopoietic structures.

## Acknowledgments

Thank to Thomas Schmickl who suggested the idea to keep track of this personal
work on arxiv.

## References

  * [1] Chou H.H., Reggia J.A. Emergence of self-replicating structures in a cellular automata space. Physica D: Nonlinear Phenomena, 110, 252-276 (1997). https://doi.org/10.1016/S0167-2789(97)00132-2.
  * [2] Sayama H. A new structurally dissolvable self-reproducing loop evolving in a simple cellular automata space. Artificial Life, 5(4), 343-65 (1999). https://doi.org/10.1162/106454699568818
  * [3] Schmickl T., Stefanec M., Crailsheim K. How a life-like system emerges from a simplistic particle motion law. Scientific Reports 6, 37969 (2016). https://doi.org/10.1038/srep37969
  * [4] Hutton T.J. Evolvable self-reproducing cells in a two-dimensional artificial chemistry. Artificial Life, 13(1), 11-30 (2007). https://doi.org/10.1162/artl.2007.13.1.11
  * [5] Luisi P.L. Autopoiesis: a review and a reappraisal. Naturwissenschaften 90, 49–59 (2003). https://doi.org/10.1007/s00114-002-0389-9

##  Appendix A Appendix

###  A.1 Program Code

Please refer to github repository:

https://github.com/ThomasCabaret/SimulationSkeleton/tree/arxiv (use the arxiv
branch), the readme contains the pieces of information to reproduce both
experiments.

###  A.2 Demo

A condensed trailer version of our particle system is available on YouTube:
https://www.youtube.com/watch?v=JmCN_4jTwf8