Meron-like topological spin defects in monolayer CrCl3

Noncollinear spin textures in low-dimensional magnetic systems have been studied for decades because of their extraordinary properties and promising applications derived from the chirality and topological nature. However, material realizations of topological spin states are still limited. Employing first-principles and Monte Carlo simulations, we propose that monolayer chromium trichloride (CrCl3) can be a promising candidate for observing the vortex/antivortex type of topological defects, so-called merons. The numbers of vortices and antivortices are found to be the same, maintaining an overall integer topological unit. By perturbing with external magnetic fields, we show the robustness of these meron pairs and reveal a rich phase space to tune the hybridization between the ferromagnetic order and meron-like defects. The signatures of topological excitations under external magnetic field also provide crucial information for experimental justifications. Our study predicts that two-dimensional magnets with weak spin-orbit coupling can be a promising family for realizing meron-like spin textures.

Unlike the easy-axis anisotropy, an easy-plane anisotropy yields a residual SO (2) symmetry that usually prohibits the formation of a spontaneous long-range magnetic order at finite temperature 1,2,4 . On the other hand, such a residual in-plane symmetry provides ample space for realizing the quasi-ordered states that evolve with various topological defects. A known example is the topological vortex/antivortex pairs with an algebraically decaying correlation, which can be described by the Berezinskii-Kosterlitz-Thouless (BKT) theory based on the 2D XY model 4,15 . However, spins of realistic materials own a threedimensional (3D) degree of freedom, making them different from the ideal XY model, in which the spin is confined within the 2D easy plane. Such an extra degree of freedom gives hope to many nontrivial topological spin states, such as skyrmions, magnetic bubbles, and merons 5,[16][17][18][19][20][21] .
These noncollinear spin textures can be induced by many mechanisms, such as antisymmetric exchange interactions (Dzyaloshinskii-Moriya (DM)), magnetic anisotropy, and magnetic dipolar couplings 5 . Particularly, isotropic easy-plane magnetized 2D materials can host magnetic vortices and antivortices. This type of topological defects were predicted to have curling magnetizations lying within the easy plane but pointing out of the plane around the core regions of vortices/antivortices 5 . In condensed matter physics, such a topological structure is referred to as a meron, which corresponds to one half of a skyrmion and is stabilized via pairs [22][23][24] . Merons have already been observed in various coupled magnetic discs and magnetic interfaces 18,[23][24][25][26][27] . For example, meron textures have been observed in thin plates of chiral-lattice magnet Co 8 Zn 9 Mn 3 and α-Fe 2 O 3 /Co heterostructure film at room temperature 25,27 . Recently, there was also a report about creating and stabilizing meron pairs in a continuous Py film, which inherits the in-plane anisotropy by local vortex imprinting from a Co disk 28 . However, besides magnetic discs and interfaces, a pristine single-atomic thin material that exhibits such topological defects is still absent. More importantly, the understanding of these topological defects, their interactions, and the profound relations with fundamental quantum phenomena, such as superfluidity and superconductivity, are yet unclear and forging ahead 29,30 .
In this work, we find that, because of the weak spin-orbit coupling (SOC), the magnetic dipolar interaction induced magnetic shape anisotropy (MSA) can overcome the magnetocrystalline anisotropy (MCA) to evince an easy-plane, isotropic magnetic polarization in ML CrCl3. By employing an anisotropic Heisenberg model with magnetic dipole-dipole (D-D) interactions and exchange interactions extracted from first-principles simulations, we predict the existence of meron-like topological defects in ML CrCl3. Moreover, beyond pairs of vortecies and antivortecies, higher-order states involved with more than two merons are also observed in our simulations, forming complex topological excitations. Finally, to guide experiments, we show that merons are robust against external fields, and the rich hybridization between merons and the FM order can be tuned via external magnetic fields.

Results
Easy-plane anisotropy in ML CrCl 3 . The easy-plane anisotropy is crucial for the formation of topological defects 5,31 . Usually, the anisotropy is described by the magnetic anisotropy energy (MAE), which characterizes the dependence of energy on the orientation of magnetization. There are two origins of MAE owing to the relativistic effect. The first is MCA that is mainly determined by the SOC. The second is the shape anisotropy related to MSA that is originated from the Breit modification of the relativistic two-electron energy. MSA is ascribed to the classical magnetic D-D interactions 32,33 : where μ 0 is the vacuum permeability, and m i , m j are the local magnetic moments with a spatial separation of r ij . Such a D-D interaction usually favors magnetization along the elongated direction of materials. For thin films, this provides an origin of inplane anisotropy 5,32 . Typically, the shape anisotropy can be added as a posterior term after the first-principles electronic calculations 33 . In most materials, the magnetic D-D interaction is small compared with the MCA interaction 34 . However, it may play an important role in weak-MCA magnets with negligible SOC. This leads us to study ML CrCl 3 , which is known for its weak SOC 33 . 2D CrCl 3 has been recently fabricated [35][36][37] , and the XY physics was discussed in a recent experiment on ML CrCl 3 38 . According to Eq. (1), the MSA energy is calculated by with the magnetization rotating from the in-plane direction (k) to the out-of-plane direction (?). The MCA energy is calculated by ? , which is obtained from first-principles density functional theory (DFT) + U calculations. The calculation setups and convergence tests are included in the Method section. The schematic plot of the FM anisotropy energy (taking the out-of-plane direction as a reference) is presented in Fig. 1(a). Because of the competition between MCA and MSA, the overall MAE of ML CrCl 3 is about −34 μeV (the negative sign means a preferred in-plane polarization). Thus, the zerotemperature ground-state magnetization in ML CrCl 3 is inplane ( Fig. 1(c)). Besides, the in-plane MAE is nearly perfectly isotropic, and the variation is <0.1 μeV. In this work, we choose the Hubbard U and Hund exchange J as 2.7 eV and 0.7 eV, respectively, which were used in published literatures and provided good agreements with available measurements 3,39-41 . Moreover, we have checked other choices of U and J within a reasonable range and confirmed that they do not qualitatively affect the preferred in-plane MAE and subsequent topological spin textures of ML CrCl 3 (see the Supplementary Note 1 and Supplementary Table 1). These results also agree with recent calculations and measurements 33,[35][36][37] .
For the purpose of comparison, we have done calculations for ML CrI 3 , which has a large SOC. Its MAE is around 300 μeV, resulting in a known Ising-like out-of-plane magnetic order 3,6,42,43 . The corresponding anisotropy energy and preferred zero-temperature out-of-plane magnetization are schematically plotted in Fig. 1(b, d). The calculated anisotropy energies and magnetic coupling constants are listed in Table 1 for these widely studied chromium trihalides, and ML CrCl 3 is the most promising structure to realize the in-plane ground-state polarization. Moreover, the relationship between the MCA energy and SOC strength is discussed in the Supplementary Note 2 to show that SOC is the main mechanism in deciding the MCA energy.
Low-temperature magnetic phase. To explore the magnetic phase at low temperature, we have employed a classical Monte Carlo (MC) simulation based on an anisotropic Heisenberg model with parameters from first-principles simulations. This Heisenberg Hamiltonian is essentially an approximation to the many-body electronic Hamiltonian in a localized basis, which only contains the spin degree of freedom 44 . This method provides good agreements of the Curie temperature (T c ) for ML CrI 3 and CrBr 3 with experiments 3,9,41,44,45 . More specifically, we add the magnetic D-D interaction (E DÀD ) from Eq. (1): where A describes the easy-axis single-ion anisotropy, λ 1;2 and J 1;2 represent the anisotropic and isotropic exchange couplings up to the next nearest neighbors, and B is the external magnetic field. The details of extracting those coefficients from first-principles simulations are presented in the Method section.
We begin with the intrinsic ML CrCl 3 without an external field. The evolution of the averaged magnetic moment with temperature for ML CrCl 3 is presented in Fig. 2(a). We find that the outof-plane magnetization (the blue-dotted line M z ) is almost completely quenched. Surprisingly, the in-plane magnetization (the red diamond M in ) with large fluctuations emerges even after the result is averaged over 25 ensembles. In contrast to CrCl 3 , as shown in Fig. 2(b), ML CrI 3 exhibits a normal Ising-like out-ofplane magnetism with negligible fluctuations, and the in-plane magnetization is always zero. A sharp phase transition is observed around 42 K for ML CrI 3 , which is close to the measured Curie temperature of 45 K 6 .
The inset of Fig. 2(a) further addresses the substantial fluctuations of the in-plane magnetization by comparing the statistically averaged in-plane magnetic moment with a singleround simulation. Importantly, for a single-round simulation, there is no obvious trend of the evolution of the in-plane polarization even between adjacent temperatures. Such a remarkable randomness strongly indicates the existence of a weak polarized state or a quasi-ordered phase, e.g., the vortex/antivortex states described by BKT physics 4,15,46 . The detailed comparison of averaged magnetization for ML CrCl 3 and CrI 3 with and without D-D interactions can be found in the Supplementary Note 3.
We have also checked the potential existence of quasi-order in ML CrCl 3 by analyzing the spin-spin correlation functions, which is an effective way to identify orders 15,47 . In Fig. 2(c), the spin correlation of ML CrCl 3 exhibits a quasi-long-range decay to zero at a low temperature (4 K) and an exponential decay to zero at a high-temperature (12 K) that is within the paramagnetic (PM) phase. These behaviors indicate that (1) there is no long-range order in ML CrCl 3 ; (2) quasi orders may exist at low temperature because of the algebraic-like decaying correlation function. On the other hand, ML CrI 3 exhibits a normal FM behavior in Fig. 2 Monolayer  Table 1 Anisotropy energies and magnetic coupling strengths of ML CrX 3 (X = I, Br, Cl).
The energy and coupling strength are in μeV and meV, respectively.
(d): the spin correlation exponentially decays to a nonzero value at 40 K because of the long-range FM order while it exponentially decays to zero at 48 K because of the PM phase. In the following, by further analyzing the real-space arrangement of local magnetic moments as obtained through MC simulations, we confirm the topological vortex/antivortex type of defects within a 3D spin space. Figure 3(a-c) present the realspace magnetic moments at different temperatures. At low temperatures ( Fig. 3(a, b)), we can identify the vortex/antivortex type of topological spin defects, which are marked by the ovals. At higher temperature (Fig. 3(c)), the spin structures disappear. Importantly, the numbers of vortex and antivortex defects are always the same, maintaining an overall integer topological unit. The chirality of topological pairs is indicated by the blue (vortex type) or red (antivortex type) part of the ovals. In addition, we have plotted the corresponding relative phase maps, which are used in experiments to identify spin textures 48,49 . In Fig. 3(d-f), the phase maps clarify the chirality of topological defects at core regions.
At low temperature (0.25 K), we find that these vortex (antivortex) type spin defects can be connected with neighbor alternative antivortex (vortex) type defects via a spin flux closure as indicated by the dashed line in Fig. 3(d). The enlarged Fig. 3(g) provides a clear micro-structure diagram of those topological defects. Interestingly, unlike the strict XY model where spin rotator is constrained within the 2D XY plane, the topological structures in ML CrCl 3 contain 3D information, and this 3D spin effect is particularly significant around the core regions of the topological structures. As seen from the side view in Fig. 3(h), the magnetic moments around cores point out of the plane and form alterable spin hills 26,28,[49][50][51] . Compared with the vortex and antivortex pairs in the XY-model, whose free energy mainly comes from the core regions 4,15 , the benefit of forming such 3D spin texture is manifest by avoiding large swirl angles around the cores of defects and, hence, lowering the system energy. Meanwhile, our simulation does not display any correlation between the magnetic moment swirling and the direction of the central spin hill. This agrees with the experimental observations that all four possible up/down combinations of paired meroncore polarities show up with nearly equal probabilities 28 .
Such topological defects have been referred to as merons, which have a ± 1=2 skyrmion charge 5,23,52 and are similar to experimental observations 18,[23][24][25][26] . This is a material analogy to a solution of the Yang-Mills theory, in which merons as described in the context of quark confinement can only exist in pairs owing to the one half of topological charge carried 22,23,52,53 . This character agrees with the equal numbers of vortex/antivortex defects observed in our simulations. Particularly, because of the pairing character of these topological spin defects, it is expected that doped free carriers could be intrinsically spin-polarized, and the transport of charged meron pairs might be dispersionless at low temperature [54][55][56] . Although those dynamic behaviors and quantitative interactions between merons are beyond the scope of this study, they would be interesting topics for future studies 57-59 .  Beyond meron pairs, our MC simulations reveal more complicated spin textures, such as the higher-order states involved more than two merons. As indicated in Fig. 3(a, d), some pairs of vortex/antivortex type merons are found to be antiparallelly aligned with each other at low temperature, forming "quadrupole-like" topological excitations 57,60 . The presence of these hierarchical excitations indicates that, although the classical MC treatment does not take into account the quantum nature of spins, it does include the correlations of magnetic moments from the Heisenberg model 9 . (see Supplementary Fig. 4 for more examples of coupled meron pairs).
As temperature increases, those higher-order excitations are firstly annihilated by thermal fluctuations as indicated by the snapshot at 3.25 K in Fig. 3(b, e). On the other hand, the relatively robust pairing between vortex/antivortex type merons still persists, although thermal energy blurs the topological structures and makes them harder to be recognized in Fig. 3(e). Finally, at the higher temperature (~12 K), the paired merons totally melt, as shown in Fig. 3(c, f). The system becomes a disordered, PM phase with an exponentially decaying spin-spin correlation, which agrees with Fig. 2(c).
There are a few important mechanisms that are crucial for understanding the simulation results. First, we clarify the role of D-D interactions in generating meron-like spin textures. It is known that D-D interactions complicate system behaviors due to its nonlocal character and weakened singularities of the longitudinal mode of spin waves [61][62][63] . For 3D spins in ML CrCl 3 , D-D interactions mainly impact the MAE and introduce the preferred in-plane ground-state polarization. This is evidenced by the similar averaged magnetic moments calculated with different truncations of D-D interactions in Fig. 2(e). In the above simulations, we truncate the D-D interaction with a cutoff radius around 31 Å, as adopted in previous studies 50,64 . When the truncation radius increases from 15 Å to 52 Å, the statistically averaged value of magnetic moment at a certain temperature (T = 3 K) does not change essentially. Moreover, we find that the residual D-D interaction beyond a 31 Å truncation contributes <5% to the total MSA energy. Thus, the overall D-D interaction energy is, at least, an order of magnitude smaller than the energy from the isotropic exchange interaction J. As a result, except at extremely low temperature, the residual energy due to the D-D interaction truncation shall not qualitatively affect the spin textures of our studied systems. It is worth mentioning that our simulation cannot give accurate asymptotical behaviors approaching zero temperature due to the exponentially increasing simulation time, but this limit does not affect the main results of this work.
The simulation-size effect also needs to be considered. As Bramwell and Holdsworth discussed, a spontaneous magnetization could show up in a finite-size XY model owing to the suppression of Goldstone modes that usually destroy ferromagnetism 63,65 . A similar behavior could happen in 3D spin cases. As shown in Fig. 2(f) at 3 K, when the simulation size is not big enough (less than~400 Å), the averaged magnetic moment of ML CrCl 3 remains nearly a constant and does not exhibit large fluctuations, resembling the typical FM order. This also agrees with previous anisotropic Heisenberg-model studies 64,66 on planar phases of 2D lattices, in which only the planar FM type curve was observed due to small simulation sizes. On the other hand, when the simulation size is significantly larger than that of meron pairs (~200 Å estimated from Fig. 3(a)), quasi orders start to show up because the system is large enough to hold meron-pairs. Particularly, with simulation size increasing, owing to the geometric features and growing density of merons, the averaged magnetization is expected to gradually decrease. This is what we observed in Fig. 2(f (a-c). The schematic spin directions and the connections between defects are drawn to clarify the vortex/ antivortex type meron and the higher-order excitations. g, i, j Are enlarged views of the specific defect pairs and the rectangular area from (a-c), respectively. h Side view of (g) with an enlarged z-axis ratio (10:1). NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18573-8 ARTICLE NATURE COMMUNICATIONS | (2020) 11:4724 | https://doi.org/10.1038/s41467-020-18573-8 | www.nature.com/naturecommunications averaged in-plane magnetic moment of ML CrCl 3 (the deep blue line with red circles) decreases with large fluctuations. In contrast, the magnetic polarization of Ising-like ML CrI 3 (the green line with light blue circles) does not depend on the system size, showing a normal FM state. Therefore, we conclude that a finitesize simulation can capture those meron pairs if the system size is significantly larger than that of merons. In other words, larger systems can hold more meron pairs but do not change the existence of the meron-like phase.
Finally, we also check the proximity effect from substrates, such as the widely used hexagonal boron nitride (h-BN). We do not find significant modifications to the meron states when ML CrCl 3 is attached to h-BN, as shown in the Supplementary Table 2, Supplementary Fig. 3, and Supplementary Note 4. This verifies the feasibility of experimental realization for our theoretical predictions.
Response to external field. It is highly motivated to further investigate the response of the topologically paired merons to an external static magnetic field, which has been widely employed to study magnetic properties of materials 6,7,35,37 . As discussed above, the polarization direction of the local magnetic moments is crucial for determining the existence of the meron phases. As a result, a magnetic field applied along different directions is expected to induce different responses. First, we apply a small inplane magnetic field (70 Gs), and the results are shown in Fig. 4 (a). At low temperature (below 2 K), the magnetic moment is close to the saturated value (3 μ B /Cr), and the fluctuation is substantially quenched, indicating a well-defined FM order. This is confirmed by Fig. 4(c, g), in which both the real-space magnetic moment distribution and projected magnetic component along the external-field direction show a FM order. This suggests that the applied in-plane magnetic field breaks the in-plane isotropy, forming a preferred polarization direction, which stabilizes the FM order. Interestingly, as temperature goes up, the magnetic polarization decreases, and larger fluctuation accompanies. In the realspace snapshot (Fig. 4(d)), meron pairs emerge. Those thermally excited topological defects are embedded in the aligned spin sea, which is tilted along the direction of the external field. Figure 4(d, h) provide a physics picture analogy of spin distributions in the real space and projected space. The latter can be directly compared with the experimentally observed figure format 48,49 .
We then explore the case where the applied magnetic field is perpendicular to the material plane. If the out-of-plane magnetic field is below a critical value (B crit ? % 0:2T), although it is not strong enough to completely govern the preferred magnetization direction, the local magnetic moments are no longer perfectly within the material plane. Consequently, we observe a net averaged magnetization along the out-of-plane direction. In this weak magnetic field situation, the topological meron pairs may be preserved. Take B ? ¼ 0:1T as an example shown in Fig. 4(b). At low temperature, an averaged out-of-plane magnetization shows up with the help of external field (the orange-dotted curve M z ). Meanwhile, the in-plane magnetic moment (the red-dotted curve M xy ) exhibits large fluctuations, indicating the existence of excited merons. This is evidenced by both the real-space and projected component plots of the magnetic moments in Fig. 4(e, i). The hot spots in the projected component figure indicate the core regions of merons, which correspond to the previously discussed spin hills at meron cores. When the temperature reaches 10 K, the averaged in-plane magnetization approaches zero. On the other hand, the out-of-plane magnetization is preserved, as shown in Fig. 4(f, j). This is another type of hybridization between FM and meron states. Moreover, the reason for the abnormal elevation of  in (a, b). The lower panels (g-j) are the corresponding projection components of local magnetic moments along external field directions (parallel and perpendicular to the material plane, respectively) from (c-f). The error bars are the standard deviation.
out-of-plane FM polarization with increasing temperature at around 10 K is that the D-D interaction is more sensitive to the arrangements of adjacent magnetic moments than the on-site Zeeman energy provided by the external field.

Discussion
In summary, our studies show that the overall isotropic in-plane magnetized 2D material such as ML CrCl 3 could provide the opportunity for generating meron-type topological defects without involving any other interactions. This would give hope to realizing topological spin structures in a much broader range of materials with weak SOC, such as CrF 3 . Meanwhile, these meron-like topological spin defects may shed light on the understanding of the recently observed suppressed magnetic ordering in ML NiPS 3 , an AFM easy-plane anisotropy material 67 . In addition to the previously studied Josephson junction arrays or superfluid helium systems. 15,29,54 , these meron-like pairs may provide an alternative way to study the exotic relations between topological spin textures and quantum critical phenomena, including the superfluidity and superconducting behaviors, in pristine 2D magnetic systems.

Methods
Density functional theory (DFT) calculations. The DFT calculations are performed within the generalized gradient approximation using the Perdew-Burke-Ernzerhof functional as implemented in the Vienna Ab initio Simulation Package 68 . A plane-wave basis set with a kinetic energy cutoff of 450 eV and a 5 × 5 × 1 k-point sampling is adopted for a 2 × 2 × 1 supercell to mimic different magnetic configurations for extracting magnetic interaction constants. The k-point grid is 12 × 12 × 1 for a 2 × 2 × 1 supercell. The vacuum distance is set to be 20 Å between adjacent layers to avoid spurious interactions. The van der Waals interaction is included by the DFT-D2 method 69 , and SOC is considered. We choose the Hubbard U = 2.7 eV and Hund J = 0.7 eV for Cr 3+ ions. The structure is relaxed until the force converges within 0.01 eV/Å. The optimal lattice constant is 6.01 Å, which is close to the bulk monoclinic structure lattice constant 6.06 Å 70 .
Determination of the Heisenberg Hamiltonian coefficients. We obtain the magnetic interaction coefficients of the anisotropic Heisenberg model Hamiltonian (Eq. (2)) by calculating the total energies of different magnetic configurations. We consider the FM and Néel antiferromagnetic (AFM) configurations. The corresponding energies for a unit cell are in which two magnetic orientations (in-plane and out-of-plane) are calculated to determine the anisotropic coupling constants. Moreover, we flip the magnetic moment of a Cr 3+ cation in the 2 × 2 × 1 supercell to realize more magnetic configurations (the energies are normalized to one unit cell).
As a result, we obtain six equations, which are used to solve for the five magnetic interaction coefficients in Eq. (1) and the reference energy (E 0 ). The extracted coefficients and corresponding magnetic anisotropy energies for CrX 3 (X = Cr, Br, I) are summarized in the Table 1.
MC simulations. Based on the Hamiltonian Eq. (2), we perform MC simulations via the Metropolis algorithm on 2D hexagonal lattices with a size of 160 × 160 unit cells (if not particularly specified), which contain 51,200 magnetic moments. The periodic boundary condition is implemented. All magnetic moments are set to along the out-of-plane direction at the initial state to mimic experimental cooling conditions with the help of external field. We run for 5:12 10 9 MC steps (10 10 5 steps per site in average) to ensure that the thermal equilibrium is achieved.
The out-of-plane and in-plane magnetizations are defined as: where N represents the number of total magnetic moments in the simulated system, and <X> denotes the time average of corresponding components after the system achieves the thermal equilibrium. In magnetic systems, the static spin-spin correlation function describes the average scalar product of spins at two lattice sites with a fixed distance r 1 À r 2 j j¼ r. It can be written as: Data availability The data that support the findings of this study are available from the corresponding author upon reasonable request.