Let us first consider the crystal structure of monolayer Cr2Ge2Te6, as illustrated in Fig. 1. The crystal adopts the space group P\(\stackrel{-}{3}1m\). Obviously, monolayer Cr2Ge2Te6 presents overall inversion symmetry, and thus the DMI between the first nearest neighbors (NN) will cancel out upon space inversion, in accordance with Moriya’s first rule [20]. However, it is noteworthy that local inversion symmetry is inherently broken for the second nearest neighbors. Let’s take Cr1-Cr2 pair as an example. As indicated in Fig. 1, there is a two-fold rotation axis C2y perpendicular to the line connecting the second NN. According to Moriya’s fourth rule [20], the DM vector is perpendicular to this axis. Consequently, the DM vector associated with Cr1-Cr2 pair possesses a z component (dz) and a planar component (dx) aligned along Cr1-Cr2 bond. Considering the three-fold rotational symmetry around [001] and the translational symmetry, the in-plane component of the DM vector between Cr1 and Cri (i = 2–7) points towards (or away from) Cr1. While the out-of-plane component alternates in sign every 60 degrees, specifically, the orientation of dz component remains consistent for Cr1-Cr2, Cr1-Cr4, and Cr1-Cr6 pairs, which is opposite to that of the Cr1-Cr3, Cr1-Cr5, and Cr1-Cr7 pairs.
Based on the symmetry analysis presented above, we demonstrate that a finite second NN DMI is allowed to exist in centrosymmetric 2D magnetic materials with a honeycomb lattice owing to the broken local inversion symmetry. Usually, the ratio of first NN DMI (d1) to the J1 is regarded as an important parameter to assess the feasibility of forming topological spin textures. The typical range for |d1/J1| falls within 0.1 − 0.2 [31]. However, for monolayer Cr2Ge2Te6 and other 2D magnetic materials with similar structure, the d1 is inhibited by inversion symmetry. To deeply understand this unusual skyrmionics physics in centrosymmetric 2D magnetic materials with a honeycomb lattice, we thus introduce the dimensionless parameter\(\kappa ={(\frac{4}{\pi })^2}\frac{{AK}}{{{D^2}}}\), that can help us identify magnetic ground state [32, 33]. Here A is the exchange stiffness that includes both Heisenberg exchange interaction J1 and J2, K is the total magnetic anisotropy, and D is the DMI parameter in micromagnetic format, respectively. It has been demonstrated that, when 0 < κ < 1, the magnetic ground state of the system exhibits spin spiral. Under a moderate magnetic field, spin spiral might transform into skyrmions. Whereas for κ > 1, the system tends to exhibit a collinear magnetic structure, i.e., a FM state.
The material parameter A is related to the parameters of the atomistic model by:
$$A=\frac{1}{{2V}}{\sum\limits_{j} {{J_{0j}}(R_{{0j}}^{x})} ^2}{\text{=}}\frac{{{a^2}}}{{2V}}(\frac{1}{2}{J_1}+3{J_2})$$
1
Here\(J_{{0j}}^{{}}\)is the NN and second NN Heisenberg exchange interactions, \(R_{{0j}}^{x}\) is the x component of Cr pair vector, respectively. V is the volume of the magnetic part of the unit cell, e.g. the Cr2Ge2Te6 monolayer, and a is the distance between two second NN Cr atoms. As indicated by our symmetry analysis, the orientation of dz for adjacent second NNs is opposite for centrosymmetric 2D magnetic materials with a honeycomb lattice, thus only the dx is left to contribute to the material parameter D. Hence, D can be expressed as follow:
$$D=\frac{1}{V}\sum\limits_{j} {d_{{0j}}^{x}R_{{0j}}^{x}} =\frac{a}{V}(3d_{2}^{x})$$
2
Finally, we can rewrite κ with the parameters of atomistic model as:
$$\kappa ={(\frac{4}{\pi })^2}\frac{{({J_1}+6{J_2})K}}{{{{(6d_{2}^{x})}^2}}}$$
3
We now examine the SIA-dx phase diagram at 0 K without external magnetic field. To this end, we adopt the following spin Hamiltonian in parallel tempering Monte Carlo (PTMC) simulations ((see Supplemental Material (SM) [34]) for details):
$$H= - \frac{1}{2}\sum\limits_{{<i,j>}} {{J_{ij}}{{\mathbf{S}}_i} \cdot {{\mathbf{S}}_j}} - \frac{1}{2}\sum\limits_{{<i,j>}} {{\mathbf{d}}_{{ij}}^{x} \cdot ({{\mathbf{S}}_i} \times {{\mathbf{S}}_j})} - \sum\limits_{i} {K{{(S_{i}^{z})}^2}} - B\sum\limits_{i} {S_{i}^{z}}$$
4
here Si is the unit vector of spin at site i. The isotropic Heisenberg exchange coupling, denoted by Jij, spans over all first and second nearest-neighbor (NN) magnetic ions. Second term \({\mathbf{d}}_{{ij}}^{x}\) represents the vector that characterizes the x component of DMI for the second NN magnetic ions. The third term is the single ion anisotropy K. In the last term, B indicates the strength of the external magnetic field along Z axis, which is same as indicting in Fig. 1. The PTMC simulations are performed over a 60 × 60 × 1 supercell. The phase diagram is obtained by fixing J1 = 10 meV and J2 = 0 meV as a function of the dx and single ion anisotropy (SIA). To obtain the ground state at 0 K, various periods of spin spirals are inputted as initial states for each set of magnetic parameters and fully relax them. The energy of the FM state (FM-x for easy-plane anisotropy, and FM-z for easy-axis anisotropy) is set as the zero reference. One can clearly see in Fig. 2 that, within the |κ| < 1 threshold, the spin spiral is indeed more energetically favorable compared to FM state for both easy-axis and easy-plane anisotropy. Even a small second NN dx, when coupled with sufficiently low SIA, results in 0 < |κ| < 1, maintaining a spin spiral ground state.
To confirm the analysis provided above, by combining DFT calculations and the energy mapping method (see SM [34] for details), we extract magnetic coefficients of monolayer Cr2Ge2Te6, including Heisenberg exchange coupling J, DMI, and magnetocrystalline anisotropy (C-MAE). The first and second NN Heisenberg exchange coupling J1 and J2 is calculated to be 10.85 and − 0.34 meV, respectively, which agree well with previous studies [35–37]. Due to the inversion symmetry, the first NN DMI is determined to be zero. In line with the aforementioned analysis, the second NN Cr pair does exhibit finite DMI. Specifically, for Cr1-Cr2 pair as depicted in Fig. 1, the planar component dx is calculated to be 0.27 meV, the perpendicular component dz is -0.61 meV, while the other planar component dy is zero. Magnetic anisotropy energy (MAE) is also a crucial factor that influences the formation of topological spin texture. The MAE is composed of magnetocrystalline anisotropy energy (C-MAE) and magnetic dipolar anisotropy energy (D-MAE). The calculated C-MAE is 0.065 meV per Cr atom (prefers an out-of-plane magnetization), and the D-MAE is -0.039 meV per Cr atom (favors an in-plane magnetization). Thus, the total magnetic anisotropy energy is out-of-plane for monolayer Cr2Ge2Te6, which is also consistent with both theoretical and experimental works [38, 39].
After obtained all the magnetic parameters, we now investigate the spin textures of monolayer Cr2Ge2Te6 under different magnetic fields at 0 K. The spin Hamiltonian of Eq. (4) is used in PTMC simulations. For simplicity, both the C-MAE and D-MAE are included into the SIA. The PTMC simulations are also performed over a 60 × 60 × 1 supercell. The initial spin configurations are obtained through fully relaxing the paramagnetic (random) state. To ensure the converged spin structures to be the ground state at 0 K, the conjugate gradient (CG) method is subsequently applied to relax the spin configurations obtained through PTMC simulations. Remarkably, without applying magnetic field, Bloch-type labyrinth domains are observed in monolayer Cr2Ge2Te6 at 0 K, as depicted in Fig. 3(a, c). For monolayer Cr2Ge2Te6, the κ is determined to be 0.288, indicating a spin spiral ground state, which aligns closely with the results from PTMC simulations. Although there are many 2D magnetic materials that share the similar structure as Cr2Ge2Te6, most of them exhibit collinear magnetic structure. Our phase diagram can explain this very well. It is known that Cr2Si2Te6 is an FM insulator. Our calculated magnetic parameters lead to a κ value of 3.94, also suggesting FM ground state (see Table SI of the SM [34]). Likewise, using the magnetic parameters available for monolayer CrI3 [40], its κ is determined to be 25.84, indicating again FM ground state.
More intriguingly, Bloch-type skyrmion is observed on the background of labyrinth domains at 0.17 T, as shown in Fig. 3(a). The sense of in-plane spin rotation for both labyrinth domain and skyrmions are exactly consistent with the experimental observations [25, 27]. When the magnetic field is larger than 0.2 T, the labyrinth domains are totally wiped out. Importantly, as illustrated in Fig. 3 (b), this skyrmion state is preserved within the range of 0.2–0.55 T. Above 0.55 T, the monolayer Cr2Ge2Te6 is completely magnetized into the trivial ferromagnetic (FM) state. It should be noted that the skyrmions state is metastable within the field range of 0.31–0.55 T, which possesses a higher energy compared to the corresponding FM state. We note that Khela et al. [26] have also proposed that the second NN DMI plays a crucial role in the formation of skyrmions in Cr2Ge2Te6. However, in their atomistic spin dynamics simulations, three distinct types of topological spin texture are found: Bloch-type skyrmions and anti-skyrmions, and skyrmioniums. Importantly, the sense of in-plane spin rotation of the skyrmion identified in their simulations however contradicts the previous experimental results, which means the orientation of dx component of second NN DMI used in their simulations is opposite to the direction derived from our DFT results.
Finally, we establish anisotropy versus magnetic field phase diagram at 0K, as a function of dimensionless magnetic field HJ/D2 and anisotropy KJ/D2, with fixed D (dx) = 0.3 and J = 10 meV, respectively. To construct the phase diagram (K, H), we compare the energy of the states obtained by PTMC to that of the uniform (tilted) FM state, with the assistance of magnetization (mz) and topological number. The uniform FM state can either be the spin polarized FM state along the magnetic field direction (nz = 1) or the spin tilted FM state with nz < 1. As shown in Fig. 4, surprisingly, the easy-plane regime (K > 0) leads to an unexpectedly large stable SkX phase, compared to the easy-axis anisotropy (K < 0). This large region of the skyrmion/skyrmions crystal (SkX) within the phase diagram predominantly aligns around dashed line H = 2A in Fig. 4, marking the boundary between the "tilted FM" and the easy-axis FM. As demonstrated in the previous study [41, 42], this is due to: (i) the rotated spins within the skyrmion configuration reduce the DMI contribution to the free energy as compared to the FM state for easy-plane anisotropy; (ii) the skyrmion presents a more optimal balance between an easy-plane anisotropy and a field along the z-axis compared to a spiral configuration. Note that our PTMC calculations fail to give clear results near the phase boundaries because the system gets trapped in metastable states. Hence, the phase diagram depicted in Fig. 4 should be regarded only as semi-quantitative.