**Spin Hamiltonian and DFT Results**

The crystal structure of monolayer CrSBr is illustrated in Fig. 1a. The lattice constants, optimized by Density Functional Theory (DFT) (see **Methods**), are *a* = 3.506 Å, *b* = 4.698 Å, respectively. Monolayer CrSBr belongs to the centrosymmetric *Pmmn* space group, the inversion center being located at the center of second-nearest (2nd, label 0 and 2 in Fig. 1a) neighbor Cr-Cr pairs. However, due to the local site inversion breaking, no inversion center occurs for the first-nearest (1st, 0 and 1) and third-nearest (3rd, 0 and 2) neighbor Cr-Cr pairs. Therefore, according to Moriya rules22, the DM interaction22,23 is along the **b** (**a**) direction and labeled as *D**y* (*D**x*) for 1st (3rd ) Cr-Cr pair in the following.

The magnetic interactions are modeled by using the spin Hamiltonian with the following general form:

$$H=\frac{1}{2}{\sum }_{i\ne j}{S}_{i}{J}_{ij}{S}_{j}+{\sum }_{i}{S}_{i}{A}_{i}{S}_{i}$$

1

where *J**ij* and *A**i* represent the exchange coupling interaction tensor and single ion anisotropy (SIA), respectively, where all the parameters can be obtained from DFT calculations using the four-state method. The *J**ij* results are shown in **Supplementary Table 1**. The diagonal terms *J**xx*, *J**yy*, and *J**zz* of the three nearest neighbors are all positive, showing ferromagnetism and preventing magnetic frustration in monolayer CrSBr. The antisymmetric part of *J**ij* corresponds to the DM interaction and can be expressed as \({D}_{\gamma }={\frac{1}{2}\epsilon }_{\alpha \beta \gamma }({J}_{\alpha \beta }-{J}_{\beta \alpha })\), where \({\epsilon }_{\alpha \beta \gamma }\) denotes the Levi-Civita symbol. The DM interactions of 1st and 3rd nearest Cr-Cr pairs in CrSBr are *D**y* = 0.16 meV and *D**x* = -0.29 meV, respectively, in agreement with the symmetry analysis by Moriya rules. Moreover, the third nearest-neighbor exchange constant, \({J}_{03}=({J}_{xx}+{J}_{yy}+{J}_{zz})/3\), is 2.96 meV and *D**x*/*J*03 is ~ 0.1, i.e., in the typical range for the formation of topological spin textures24. In addition, due to the orthorhombic structure of CrSBr, the *A**xx*, *A**yy*, and *A**zz* parameters are not equal to each other; rather, the calculated *A**xx*-*A**zz* and *A**yy*-*A**zz* are − 0.01 meV and − 0.02 meV, respectively, consistent with experiments suggesting the easy, hard, and intermediate axis of CrSBr to be along the **b**, **c**, and **a** direction, respectively 15.

**Equilibrium Monte-Carlo Simulations**

Based on Eq. (1), we performed atomistic spin simulations via Monte-Carlo (MC) (see **Methods**) for a 100 nm × 100 nm monolayer CrSBr, focusing on its magnetic transition temperature. On one hand, a classic FM *M-T* curve is shown in Fig. 1b, indicating a FM ground state of monolayer CrSBr with Curie temperature *T**c* ~ 112 K. The *T**c* is slightly lower than experiments, possibly due to computational parameters (such as the Hubbard U value within DFT + U) or to the underestimated lattice constant25. To determine the magnetic order universality class, we fit the curve using the relation \(M={M}_{0}{[1-\frac{T}{{T}_{c}}]}^{\beta }\), where *β* is the critical exponent26,27. The 2D-Ising and 2D-XY models are also plotted in Fig. 1b for comparison. The fitted *β* is 0.22, suggesting that, although the monolayer CrSBr shows a (weakly) uniaxial in-plane magnetic anisotropy, its magnetic order closely follows the 2D-XY model. On the other hand, the experimentally observed phase transition at *T** around 30 ~ 40 K is not observed in our MC simulations. This suggests that *T** might be related to some metastable states, not detected by means of equilibrium MC simulations.

**Non-equilibrium Spin Dynamics Simulations**

In order to probe the possible emergence of an unconventional order in monolayer CrSBr at *T* < *T**, we perform non-equilibrium spin dynamics simulations via the Landau-Lifshitz-Gilbert (LLG) equation, with 10-times averaged results shown in Fig. 1c. The in-plane magnetization, \({M}_{xy}=<\sqrt{{\left({\sum }_{i=1}^{N}{S}_{i}^{x}\right)}^{2}+{\left({\sum }_{i=1}^{N}{S}_{i}^{y}\right)}^{2}}>/N\), is much higher than the out-of-plane component \({M}_{z}=<{\sum }_{i=1}^{N}{S}_{i}^{z}>/N\), indicating the easy-plane magnetism of CrSBr. However, the results of non-equilibrium spin dynamics are qualitatively different from equilibrium MC simulations, especially for the behavior at low *T* where *M/M**s* is far from 1. To have some further insights, the spin-spin correlation function, \(C\left(r\right)=<{S}_{1}\left({r}_{1}\right)\bullet {S}_{2}\left({ r}_{2}\right)>/{\left|S\right|}^{2}\), at different temperatures is evaluated by both MC (Fig. 1d) and spin dynamics (Fig. 1e) methods. For MC calculations, the correlation function sharply decays from 1 at short-range to a constant at long-range for different temperatures, with increasing magnitude upon decreasing *T*. For spin dynamics simulations, the behavior of the correlation function at high *T* (100 and 150 K) is similar to that obtained via MC. However, at intermediate *T* (40 and 80 K), the correlation function exponentially decays to nonzero values. Interestingly, at low *T* (0, 5, 20 K and *T* < 28 K in **Supplementary Fig. 1**), an algebraic-like decaying correlation function is observed, suggestive of a BKT behavior where “quasi-order” dominated by vortex-antivortex pairs may arise. We will focus on the spin dynamics results in the following, while the comparison between equilibrium MC and non-equilibrium spin dynamics calculations will be discussed at length below.

To gain a more comprehensive understanding of the magnetism at different temperatures, we further investigate the spin textures in real-space obtained with spin dynamics. The in-plane and out-of-plane magnetization are shown in the upper and lower panels of Fig. 2a, respectively. It is clear that, upon cooling, the long-range random disordered state at 150 K (see also zoom in Fig. 2b) changes to a nearly short-range ordered magnetic-domains state at 50 K; subsequently, some magnetic defects emerge at 10 K. Some of the defects annihilate during the cooling process and, at last, some clear magnetic defects appear to be randomly distributed at domain wall at very low temperatures (the detailed spin evolution process is shown in **Supplementary Fig. 2**). The representative defects in the green and orange circular regions are shown enlarged in Fig. 2c, where in total four different topological spin textures are observed. Unlike the original 2D-XY model, where all the spins are constrained to be in-plane, our topological spin textures exhibit a 3D distribution, with spins at the core pointing out-of-plane and spins at the boundary lying in-plane. Consistent with what is expected for in-plane anisotropy, such topological defects are merons or antimerons, each having a ± 1/2 topological number28. Interestingly, due to the conserved winding number in in-plane magnets, the core-up (core-down) antivortex and core-down (core-up) vortex always occur in pair29. As a result, topological defect pairs form “bimerons”30 in monolayer CrSBr. We note here that after a long-time evolution, the bimerons disappear and the phase collapses to a trivial FM ground state. The metastable magnetic defect state is estimated to be of the order of only 0.15 meV/nm2 higher in energy than the FM state. Moreover, within spin dynamics simulations, the lifetime of the defects is ~ 1.5 ns, i.e., of the same order of magnitude as other magnetic defects experimentally observed31–33. The bimerons also occasionally emerge when using MC methods (see **Supplementary Fig. 3**).

**Influence of Different Exchange Interactions**

In order to achieve long-range magnetic order, anisotropic magnetic interactions must be included to overcome the Mermin-Wagner theorem. To further investigate the role of the different anisotropic magnetic interactions, we artificially change the value of some interactions, with respect to the DFT-derived values considered so far. First, artificial magnetic parameters that are similar to those calculated for CrSBr (i.e., *D**x* = -0.4 meV, *D**y*/*D**x* = -0.5, *J**1* = *J**2* = *J**3* = *J*, *J*/*D**x* = -6) are used to clarify the role of SIA under an out-of-plane external magnetic field *B**z*, as shown in Fig. 3. For magnets with easy-plane anisotropy (SIA < 0), the bimerons are present under low *B**z*. It is noteworthy that these bimerons exist even at SIA = 0 (see Fig. 3a and **c**). When simply setting the SIA = 0 and keep other parameters as derived from DFT, the bimerons in the domain wall region are also observed (see **Supplementary Fig. 4**). Under a suitably applied external *B**z*, the spins will tend to rotate toward the *B**z* direction, destroying the bimerons and forming the bubbles (cf. Figure 3b**)**. When *B**z* > 1T, all the spins are along the *z* direction and a trivial FM state occurs. Interestingly, as the conditions change from SIA < 0 to SIA > 0 (*i.e.*, easy-axis anisotropy), some domains occur at a smaller *B**z* (see Fig. 3). Furthermore, a suitable *B**z* (~ 0.1 T) can change a “trivial-domain” state to a “nontrivial skyrmion bubble”34 state (see Fig. 3d). As expected, all spins rotate to the *z* direction when *B**z* is strong enough (here *B**z* > 0.2 T). Therefore, within the phase diagram shown in Fig. 3, we can conclude that the bimerons emerge at SIA ≤ 0 at low *B*z. We also remark that the dipole-dipole interaction favors the in-plane spin direction in magnetic films and is reported to contribute to topological defects in CrCl335,36. We note that the tiny dipole-dipole interaction in CrSBr may increase the in-plane magnetic anisotropy; however, our results show that the bimerons always occur, no matter whether the dipole-dipole interaction is explicitly considered or not (see **Supplementary Fig. 5**). Therefore, up to this point, neither the SIA nor the dipole-dipole interaction represent the main reason for driving the bimerons formation.

Let’s now turn to another anisotropic magnetic interaction in CrSBr, *i.e.*, the antisymmetric DM interaction. Interestingly, as we artificially remove the DM interactions in CrSBr, the topological bimerons disappear, only some domains are present (see **Supplementary Fig. 6**). To further illustrate the role of DM interaction, we fix *D**x* = -0.4 meV, *J**1* = *J**2* = *J**3* = *J*, SIA = 0 and focus on the variation of bimeron states (meron-antimeron pairs) as a function of |*D**y*/*D**x*| and |*J*/*D**x*|, as shown in Fig. 4. The bimerons emerge in all conditions in this phase diagram. Furthermore, the DM interaction will favor orthogonal spins, therefore resulting in an increase of the density of bimerons, as |*J/D**x*| (|*D**y**/D**x*|) decreases (increases). Moreover, not only the number but also the lifetime of bimerons is increased for a higher *D**x* and *D**y*. For example, the lifetime of bimerons in the red circle region in Fig. 4 is ~ 75 ns and the metastable bimeron state is 0.13 meV/nm2 higher in energy than the FM state. These results demonstrate that the antisymmetric DM interaction is crucial for the formation of topological bimerons defects in monolayer CrSBr.

**Kibble-Zurek Mechanism**

As shown in **Supplementary Fig. 2**, the bimerons occur at the domain walls and annihilate one by one if, after cooling down the temperature, the system is allowed to evolve for a long time. Interestingly, as shown in Fig. 5a, the density of bimerons becomes higher as the cooling rate increases. The nearly power-law relation between the density of topological defects and the cooling rate in monolayer CrSBr is suggestive of the Kibble-Zurek (KZ) mechanism11,12. Though the KZ mechanism was first predicted in the early Universe and was later tested in liquid crystals, superfluid helium, superconductors, ultracold gases, etc., its investigation in magnets is still limited37–39. The KZ mechanism allows for the formation of defects near the phase transition; however, the cooling rate has a large influence on the evolution of the defects, suggesting that the presence of bimerons might depend on the “history” of the sample. *T** might be even changed under different cooling rates. Although we are not able to precisely determine the numerical *T** value, as the cooling rate becomes slower, there seems to be a tendency for the transition from an exponential-like to an algebraic-like decaying curve for spin-spin correlations to occur at a lower temperature (see **Supplementary Fig. 7**).

Finally, finite-size effects are considered. Figure 5b illustrates that at ~ 0 K, when the simulation size is smaller than 30 nm, both the in-plane and out-of-plane magnetizations, *M**xy* and *M**z*, show a nearly constant value and without large fluctuations, even after averaging the results 10 times. The *M**xy* ~ 1 indicates monolayer CrSBr to resemble an in-plane ferromagnet for a small system size. When the system size reaches ~ 50 nm, the magnetization drops and, importantly, large magnetic fluctuations occur. This is attributed to the size of the system being close to the bimeron size (~ 38 nm, see **Supplementary Note 1**), and both the bimerons and domains have the possibility to emerge when considering different initial random states (**Supplementary Fig. 8**). At last, when the system size is large enough, the bimerons emerge and can be uniformly distributed at the domain wall regions, leading both the averaged *M**xy* and *M**z* values to gradually decrease as the system size increases. Meanwhile, due to the KZ influence, the bimeron size changes to ~ 28nm and ~ 60nm when the cooling rate becomes two times faster and slower, respectively (see **Supplementary Figs. 9 and 10**).