Skyrmions in twisted van der Waals magnets

The discovery of ferromagnetism in two-dimensional (2D) monolayer CrI₃ and other 2D van der Waals (vdW) materials opened a new window for exploring low dimensional magnetism and its applications in spintronicsGibertini et al. (2019); Gong and Zhang (2019); Burch et al. (2018); Huang et al. (2017); Gong et al. (2017); Deng et al. (2018); Bonilla et al. (2018); O’Hara et al. (2018); Soriano et al. (2019). The properties of 2D materials can be controlled by external parametersCristoloveanu et al. (2019); Ye et al. (2012); Cao et al. (2018a); Huang et al. (2018); Song et al. (2019) and are highly sensitive to stacking and twisting between the layersPonomarenko et al. (2013); Dean et al. (2013); Gorbachev et al. (2014); Cao et al. (2018b, a). In particular, with the discovery of superconductivity in twisted bilayer grapheneCao et al. (2018b, a), there has been tremendous progress on exploring moiré superlattices both experimentally and theoretically Zhang et al. (2017); Moon and Koshino (2014); Chen et al. (2019a, b); Wallbank et al. (2015); Ni et al. (2015). In terms of magnetism, stacking order and twisting can significantly alter the interlayer exchange as the exchange is highly sensitive to atomic registriesChen et al. (2019c); Sivadas et al. (2018); Jiang et al. (2019); Du et al. (2016); Klein et al. (2018); Tong et al. (2018).

Magnetic skyrmions Roessler et al. (2006) are nanoscale vortex-like spin textures that were first observed in non-cetrosymmetic bulk magnetic materials such as MnSi Mühlbauer et al. (2009); Yu et al. (2010), (FeCo)Si Münzer et al. (2010) and FeGe Yu et al. (2011). Skyrmions are topological defects, and they exhibit novel transport phenomena such as topological Hall effect and topological Nerst effect Hamamoto et al. (2015); Shiomi et al. (2013). In recent years skyrmions received ample attention due to their potential for spintronics applications and memory storage devices Fert et al. (2013). In most cases, skyrmions are stabilized by interplay of Dzyaloshinskii-Moriya (DM) interaction and external magnetic fieldMühlbauer et al. (2009); Yu et al. (2010); Banerjee et al. (2014). In this article, we explore the possibility of stabilizing magnetic skyrmions in the absence of an external magnetic field in moiré superlattices. We consider a ferromagnetic (FM) monolayer on an antiferromagnetic (AFM) substrate with Néel order. Twisting the FM layer by an angle $\theta$ produces moiré patterns as shown in Fig. 1(a), (c). Ferromagnetic coupling between the substrate and the FM monolayer leads to an alternating exchange field for the moiré superlattice as shown in Fig. 1(b). Our setup is motivated from Ref. Tong et al., 2018. However unlike Ref. Tong et al., 2018 which includes dipole-dipole interaction to stabilize magnetic skyrmions, we consider DM interaction which is the primary interaction for magnetic skyrmions in chiral magnetsMühlbauer et al. (2009); Yu et al. (2010); Johnson et al. (1996); Dubowik (1996); Banerjee et al. (2013).

Our main results are summarized in Fig. 2 and 3. We show that (i) skyrmion crystal (SkX) is stabilized as a function of exchange coupling between the layers ($J_{\rm ex}$) and moiré periodicity. (ii) Even though SkX can be stabilized for a wide range of twisting angle, we find the optimal moiré periodicity to be about $L=9L_{D}$, where $L_{D}=(J/D)a$ is the intrinsic length scale for skyrmions. (iii) Unlike chiral magnets under magnetic field, we find an extended region of SkX for easy-axis anisotropy. (iv) We show that a large fraction of the topological charge $q=\frac{1}{4\pi}\int d^{2}r\widehat{\textbf{m}}.(\partial_{x}\widehat{\textbf{m}}\times\partial_{y}\widehat{\textbf{m}})$ of the magnetic skyrmions is concentrated at the edges and splits into three parts for large moiré periodicity and large easy axis anisotropy. This effect arises due to the anisotropic shape of the skyrmion.

Before we delve into the analysis of the effective magnetic Hamiltonian, we first describe our setup. As mentioned above, we follow the procedure of Ref. Tong et al., 2018 to derive the effective interlayer exchange field. We consider a FM monolayer twisted on top of an AFM substrate, both on a honeycomb lattice with the same lattice constant, $a$. For twisting angle $\theta$, the moiré period is given by $L=a/2\sin(\theta/2)$. For small angle $\theta$ and/or lattice mismatch $\delta$, $L\thickapprox a/\sqrt{\theta^{2}+\delta^{2}}$ (large period), the local atomic registries on length scale smaller than $L$ but larger than $a$ matches the atomic stacking of different interlayer translation r as shown in Fig. 1(c). Hence moiré superlattice can be described by interlayer translation vector r(R) that gives atomic registry at position R. The interlayer exchange coupling between AFM substrate and FM layer is different at different positions due to different atomic stacking of monolayer and substrate and this leads to spatially dependent exchange field B(R) as shown in Fig. 1(b). For example, at position $\textbf{R}_{1}$ the coupling aligns the spins of 2D layer (green) in positive z-direction when it sits on top of AFM sublattice with spins up (black) and the spins align in opposite direction at $\textbf{R}_{2}$ when it is on top of AFM sublattice with spins down (red). The interlayer exchange field at interlayer translation r is give by Tong et al. (2018)

$\displaystyle\textbf{B}^{\tau}(\textbf{r})=\sum_{j,\tau^{\prime}}J^{\tau^{\prime}}_{ex}(\textbf{r}-\bm{\xi}_{\tau}+\textbf{R}_{j})\textbf{m}_{j,\tau^{\prime}}$

(1)

where $J_{ex}$ is the interlayer coupling coefficient, m is the magnetic moment of top most layer of AFM substrate, $\left\{\tau,\tau^{\prime}\right\}=\left\{A,B\right\}$ represents the two in-equivalent sites in unit cell, $\bm{\xi}_{A}=0$, $\bm{\xi}_{B}=\left\{0,a\right\}$ and the summation $j$ is over the Bravais lattice. The total interlayer exchange field per unit cell is given by summing the fields of site $A$ and $B$

$\displaystyle\textbf{B}(\textbf{r})=\textbf{B}^{A}(\textbf{r})+\textbf{B}^{B}(\textbf{r}).$

(2)

This approximation holds when the interlayer coupling, $J_{ex}$ is small as compared to intralayer coupling, $J$. In general, $J_{ex}<J$ holds for vdW magnets as the interplane exchange is expected to be much smaller than intraplane exchange. We used the following coupling form that decays exponentially at long distances

$\displaystyle J^{\tau^{\prime}}_{ex}(\textbf{r}-\bm{\xi}_{\tau}+\textbf{R}_{j})=J^{0}_{ex}e^{-\frac{\sqrt{(\textbf{r}-\bm{\xi}_{\tau}+\bm{\xi}_{\tau^{\prime}}+\textbf{R}_{j})^{2}+d^{2}}}{r_{0}}}$

(3)

where $d$ is the interlayer separation and $r_{0}$ is the decay length. In our calculations we used $d=a$, $r_{0}=a$ and this leads to

$\displaystyle B_{Max}\thickapprox 1.9244\times 10^{-2}J^{0}_{ex}.$

(4)

Next we describe our model for the monolayer. We consider a magnetic model for 2D honeycomb lattice which is relevant to 2D vdW magnets such as trihallides Huang et al. (2017)

	$\displaystyle H$	$\displaystyle=$	$\displaystyle-J\sum_{{r,\mu}}{\bf S}_{r}.({\bf S}_{r+\hat{\delta}_{\mu}})-D\sum_{{r,\mu}}[\hat{d}_{\mu}.({\bf S}_{r}\times{\bf S}_{r+\hat{\delta}_{\mu}})]$		(5)
			$\displaystyle-A_{c}\sum_{{r,\mu}}[({\bf S}_{r}.\hat{d}_{\mu})({\bf S}_{r+\hat{\delta}_{\mu}}.\hat{d}_{\mu})]+A_{s}\sum_{{r}}(S_{r}^{z})^{2}$
			$\displaystyle-\sum_{r}{\bf B}({\bf R}_{r})\cdot{\bf S}_{r}$

where $\vec{S}_{r}$ is local moment at site $r$ and $\hat{\delta}_{\mu}$ are the three nearest neighbors on the honeycomb lattice. $J$ is the ferromagnetic Heisenberg exchange coupling, $D$ is the DM coupling Moriya (1960b). $A_{c}$ and $A_{s}$ are the compass and single-ion anisotropies respectively. The DM vector $\hat{d}_{i}=\hat{z}\times\hat{\delta}_{i}$ is set by the symmetry and originates due to the inversion symmetry breaking on the surface. However, the microscopic mechanism for the DM interaction in insulating vdW magnets is different than metallic multilayers since the former is due to the superexchange interaction in the presence of an electric field as originally discussed by MoriyaMoriya (1960b) whereas the latter can be ascribed to asymmetric interaction paths proposed by Fert and LevyLevy and Fert (1981); Fert and Levy (1980). $\bf B(R)$ is the interlayer exchange field with the twisted substrate. To explore the phase diagram of $H$, we consider the free energy functional in the continuum $\mathscr{F}[{\bf m}]=\int d^{2}rF({\bf m})$ where ${\bf m}({\bf r})$ is the local magnetization. We set the lattice constant, $a=1$. The free energy density has the following four components

$\displaystyle F({\bf m})=F_{iso}+F_{DM}+F_{aniso}+F_{moire}$

(6)

where

$\displaystyle F_{iso}$

$\displaystyle=$

$\displaystyle F_{0}(\textbf{m})+\frac{3}{2}(J/2)\sum_{{\alpha}}(\nabla m^{\alpha})^{2}$

(7)

	$\displaystyle F_{DM}$	$\displaystyle=$	$\displaystyle-\frac{3}{2}D(m^{z}\partial_{x}m^{x}-m^{x}\partial_{x}m^{z})$		(8)
			$\displaystyle+\frac{3}{2}D(m^{y}\partial_{y}m^{z}-m^{z}\partial_{y}m^{y})$

	$\displaystyle F_{aniso}$	$\displaystyle=$	$\displaystyle-\frac{3}{2}A_{c}[(m^{x})^{2}+(m^{y})^{2}]+A_{s}(m^{z})^{2}$		(9)
			$\displaystyle+\frac{3A_{c}}{4}[m^{y}(\partial_{x}m^{y}+\partial_{y}m^{x})$
			$\displaystyle-m^{x}(\partial_{x}m^{x}-\partial_{y}m^{y})]$
	$\displaystyle F_{moire}$	$\displaystyle=$	$\displaystyle-{\bf B(R)}\cdot{\bf m}.$		(10)

We absorb factor $\frac{3}{2}$ in $J$, $D$ and $A_{c}$ and define the effective anisotropy $A=A_{c}+A_{s}$ which can be positive (easy plane) or negative (easy axis or Ising). The last two terms in $F_{aniso}$, given as $\sim[m^{y}(\partial_{x}m^{y}+\partial_{y}m^{x})-m^{x}(\partial_{x}m^{x}-\partial_{y}m^{y})]$ have no contribution to the free energy for systems with periodic boundary conditions since they are total derivatives. To obtain the ground state spin configuration m, we solve the coupled Landau-Lifshitz-Gilbert (LLG) equations Gilbert (2004)

$\displaystyle\frac{d\textbf{m}}{dt}=-\gamma\textbf{m}\times\textbf{B}^{eff}+\alpha\textbf{m}\times\frac{d\textbf{m}}{dt},$

(11)

where $\textbf{B}^{eff}=-\delta H/\delta\textbf{m}$, $\gamma$ is gyromagnetic ratio and $\alpha$ is Gilbert damping coefficient. We start from different initial states¹¹1We consider multiple random configurations as well as ferromagnetic states pointing along (0,0,1), (1,1,5), (1,1,4), (1,1,3), (1,1,2), (1,1,1), (1,1,0), (2,2,1), (3,3,1), (4,4,1), (-1,-1,0), (-1,-1,-1), (-1,-1,-2), (-1,-1,-3), (0,0,-1) directions. and compare the energies of final states to get the actual ground state. To solve LLG equations numerically we used mid point method d’Aquino et al. (2005) by discretizing the effective magnetic field $\textbf{B}^{eff}=-\delta H/\delta\textbf{m}$ on a $1\times 2$ moiré supercells. We used $L_{D}=(J/D)a=10a$ to construct the phase diagrams and the results were also verified at various points for larger values of $L_{D}$. The magnitude of the magnetization was kept constant at each grid point after each time step enforcing the hard spin constraint, $|{\bf m}|^{2}=1$ and periodic boundary conditions were imposed at the boundaries. There is unavoidable error due to solving the LLG equations on a discrete lattice. The error is small as long as the variation of the magnetization is small within neighboring sites. We checked that our results are robust as a function of increasing system.

Our method can only capture magnetic orders that are commensurate with the 1x2 moiré supercell. However, in the limit of $J_{ex}\rightarrow 0$, the ground state is a spiral with the wave vector $q\simeq D/J$ for weak anisotropyBanerjee et al. (2014) which is in general not commensurate with the moiré supercell. Therefore, our method is not suitable to explore the weak $J_{ex}$ limit. Such incommensurate phases and incommensurate to commensurate phase transitions have recently been analyzed in Ref. 47. The parameters used in this article lie outside the regime of incommensurate phasesHejazi et al. (2020), yet it is not possible to extrapolate our results to $J_{ex}\rightarrow 0$ limit due to this reason.

We start by exploring the interplay between the moiré periodicity $L$ and $B_{Max}$ as shown in Fig. 2 for $A=0$. $B_{Max}$ is the maximum value of interlayer exchange magnetic field. As shown in Fig 1(b), the moiré supercell splits into two triangles of opposite alternating effective magnetic field which is maximum at the centers and vanishingly small at the corners of the triangles. Therefore, in the limit of large $B_{Max}$, the magnetization aligns with the effective magnetic field at the center of the triangles, creating two ferromagnetic domains separated with a domain wall whose chirality is set by D. The only degree of freedom left that determines the ground state is the magnetization at the corners of the triangles which creates an effective triangular lattice. In particular, if the magnetization at the corners are ferromagnetic, the state is a skyrmion whereas a stripe order give rise to a spiral phase. We find that the ground state is a spiral for low exchange field. As we increase the field, the SkX phase starts at the moiré period $L/L_{D}\approx 9$ which corresponds to the optimum angle between the layer and the substrate. As we further increase the field, we get SkX for a range of values around optimum period. This range increases with increasing the exchange field. Unlike skyrmions in chiral magnets where the size of the skyrmion is set by $L_{D}$, here we find that the size of skyrmions is determined by the moiré period. On the other hand, $L_{D}$ determines the boundary length between the interior and the exterior of the skyrmions. For small moiré period, skyrmions are small and their shape is nearly circular as shown in Fig. 2(b). Fig .2(c) shows that as the period increases, the size of skyrmion also increases and it takes the triangular shape of the exchange field. The corners of skyrmions get sharper with increasing $L$. Unlike the skyrmions in chiral magnets, we find that a large fraction of the topological charge is concentrated at the edges of the skyrmion. This fraction increases with increasing $L$. There is also a small fraction of opposite charge between the skyrmions which decreases with increasing $L$. This charge arises due to the anti-vortices between the skyrmionsLin et al. (2015). For large $L$, the topological charge further splits into three parts due to the triangular shape of the skyrmion as shown in Fig. 2(g).

Next, we explore the effects of anisotropy, $A$ at around optimal angle $L=8L_{D}$ as shown in Fig. 3. At low exchange field, we obtain spiral phase for a wide range of $A$, a small ferromagnetic phase at lowest negative values of $A$ as well as a mixed (FM+Sp) phase at largest positive values of $A$. As the field increases, initially we get SkX near the lowest negative values of $A$ that corresponds to easy axis anisotropy. By further increasing the field, the range of SkX gradually increases and eventually, occupies the whole phase diagram. Fig. 3(b,c,d) show the local magnetization and Fig. 3(e,f,g) show the spin structure factor $I(Q)\propto|\langle\textbf{m}_{Q}\rangle|^{2}$ for the three phases. $\textbf{m}_{Q}$ is the Fourier transform of the magnetization. Unlike an isotropic SkX which has six peaks on the circle in the spin structure factor, we find four peaks lie on circle and two peaks lie inside the circle. This is due to the anisotropic triangular shape of the SkX. The spiral has two peaks at $Q=\pm Q_{0}$ and the mixed state (FM+Sp) has three peaks including the $Q=0$ from the FM and $Q=\pm Q_{0}$ from the spiral phase. The intensity of $Q=0$ peak increases with increasing $A$ and decreases with increasing exchange field $B_{Max}$. On the other hand, the intensities corresponding to the spiral wave vectors have the opposite behavior of $Q=0$ with $A$ and $B_{Max}$.
The properties of SkX also depends on the anisotropy $A$. Fig. 4 shows the magnetization, spin structure factor and topological charge density as a function of $A$ for moiré period $L=8L_{D}$. For $A=-1.2D^{2}/J$ the skyrmion has a sharp boundary wall where magnetization changes abruptly and then it changes slowly inside the skyrmion. The sharpness of boundary wall decreases with increasing $A$ and the change in magnetization inside the skyrmion increases with increasing $A$. For $A=-1.2D^{2}/J$, a large fraction of topological charge is concentrated at the boundary wall of skyrmion and a small fraction lies inside the skyrmion. There is also a small fraction of opposite charge between the skyrmions. The concentration of charge at boundary decreases with increasing $A$ and the central charge increases with increasing $A$. The fraction of opposite charge between the skyrmions also increases with increasing $A$Lin et al. (2015).

We also studied the effects of anisotropy for a non-optimal angle at $L=5L_{D}$. As shown in Fig. 5, SkX is highly suppressed in this case but still persists for large exchange field and easy axis anisotropy. Suppression of SkX is due to the fact that the moiré supercell is too small with respect to the optimal size of the skyrmions.

Our results apply to a wide range of vdW magnets since the primary constraint in our model is the requirement that $J_{ex}<J$ which is in general satisfied for vdW magnets. Our phase diagrams span a wide range of parameters and for a fixed twisting angle, our results scale with $AJ/D^{2}$.

We have shown that skyrmion crystal can be stabilized in moiré superlattices in the absence of external magnetic field. We found a large SkX phase for easy axis anisotropy which can be essential to stabilize skyrmions in vdW magnets such as CrI₃Huang et al. (2017). In particular, for optimal moiré periodicity SkX occupies the majority of the phase diagram. We find that the properties of the skyrmion can be tuned with the moiré periodicity and anisotropy. Unlike skyrmions in chiral magnets, the topological charge density depends on the size of the skyrmions sand it is concentrated at the edges for skyrmions with large L.

We thank Sumilan Banerjee for useful discussions. MA is supported by Fulbright scholarship and ASU startup grant and OE is supported by NSF-DMR-1904716. We acknowledge the ASU Research Computing Center for HPC resources.