Spontaneous skyrmionic lattice from anisotropic symmetric exchange in a Ni-halide monolayer

Topological spin structures, such as magnetic skyrmions, hold great promises for data storage applications, thanks to their inherent stability. In most cases, skyrmions are stabilized by magnetic fields in non-centrosymmetric systems displaying the chiral Dzyaloshinskii-Moriya exchange interaction, while spontaneous skyrmion lattices have been reported in centrosymmetric itinerant magnets with long-range interactions. Here, a spontaneous anti-biskyrmion lattice with unique topology and chirality is predicted in the monolayer of a semiconducting and centrosymmetric metal halide, NiI2. Our first-principles and Monte Carlo simulations reveal that the anisotropies of the short-range symmetric exchange, when combined with magnetic frustration, can lead to an emergent chiral interaction that is responsible for the predicted topological spin structures. The proposed mechanism finds a prototypical manifestation in two-dimensional magnets, thus broadening the class of materials that can host spontaneous skyrmionic states.

The tensor (S1) can be diagonalized analytically, with eigenvalues (λ) and eigenvectors (ν) given by  (Table I) were calculated for the Ni 0 -Ni 1 pair, whose bonding vector is parallel to the cartesian axis x, as shown in the central panel; corresponding components for the other symmetry-equivalent pairs in the chosen cartesian reference system can be deduced via the three-fold rotoinversion symmetry associated to the D 3d point group of the analysed NiX 2 monolayers. The cartesian components of the exchange tensor for Ni 0 -Ni 3 and Ni 0 -Ni 5 bonds are therefore obtained via a rotation of the local cartesian reference system (with axis x parallel to the chosen bond) by ±120 • around the c axis, opposite pairs being related by inversion symmetry. The same symmetry arguments also apply for the second-and third-nearest neighbour interactions [shown in Fig. S2(c-e)].
canting of the two-site anisotropy axes from the perpendicular z direction. This can be even clearer by considering two extreme cases, that is either J yz = 0 or (J yy − J zz ) = 0 -i.e. isotropic diagonal terms J xx = J yy = J zz = J -in the above expressions : implying coplanar principal axes on the Ni triangular lattice and locally parallel to the x, y, z cartesian axes, i.e. spins interaction independent on the Ni-I-Ni-I plaquettes, as shown in the inset (ν α -red, ν β -green, ν γ -blue). In this case, there are no terms able to drive the formation of a non-coplanar spin-texture, which can then result into a trivial non-collinear spin-configuration, still driven by the magnetic frustration, as it results from performed MC simulations test.
which still introduces non-coplanarity of principal axis, as shown in the inset; nevertheless, having reduced the strenght of the global anisotropy in terms of energy, the formation of the non-coplanar spin-texture, i.e. the A2Sk, would depend on the competition between the isotropic term and the anisotropic J yz . In fact, from MC simulations test (not reported here) in this extreme case an increased value of J yz is required to obtain the spontaneous formation of the A2Sk-lattice.
Moreover, the sign of the off-diagonal J yz component, not affecting the principal values that depend on J 2 yz , determines the direction of the principal axes, being responsible for the change of sign of eigenvectors components (S3) and, consequently, for the change of chirality discussed in the main text and showed in Fig. S4.
In the basis of the principal axes, the exchange tensor for a given spin pair is diagonal: The nearest-neighbour exchange-coupling Hamiltonian can be written in such basis as: where the local orthogonal basis {αβγ} depends on each spin pair; being the principal axes basis orthogonal and assuming λ α = λ β (in agreement with the estimated values reported in Table I of the main text, and neglecting the small difference between λ α and λ β found for NiI 2 ), the Hamiltonian can be equivalently expressed as where J = (λ α + λ β )/2 can be seen as an isotropic exchange parameter and K = (λ γ − J ) parametrizes the Kitaevlike anisotropic exchange interaction coefficient, as in Refs [1,2]. Using the estimated values given in Table I of the main text, numerical values are J = −5.10, −6.0, −8.1 meV and K = 0.0, 0.3, 3.3 meV for NiCl 2 , NiBr 2 and NiI 2 monolayers, respectively. Such decomposition further confirms that the anisotropic contribution is not negligible in NiBr 2 and NiI 2 , being far more pronounced in the latter system.  (Table I), keeping the negative sign of the Jyz terms in the Ni 0 -Ni 1 pair (a) and changing its sign (b). Change in the rotational sense of spins, i.e. the chirality of the anti-biskyrmion, is observed. (c,d) In-plane projection for each Ni 0 -Ni i pair of the noncoplanar να eigenvector from the original calculated first-neighbor exchange interaction with negative Jyz (c), and from the first-neighbor exchange tensor where Jyz is imposed artificially to be positive (d). The applied change produces inversion in the (x, y) components of the eigenvector, that corresponds to reflection with respect to the {x,y} plane of an axial vector. Spins on the nearest-neighbour Ni atoms surrounding the central magnetic site of the topological core orient following the magnetic interaction, producing then the change of spins rotational sense also for the next-nearest neighbours, as highlighted in the zoomed view of the obtained spin configuration. A change of the chirality tuned by the sign of the Jyz term is also observed in the skyrmion lattice spin-texture obtained under an external magnetic field (Bz/J 1iso 1.3) (e,f ): helicity changes from η = −π/2 (e) to η = π/2 (f).

FIG. S5:
Dependence of vorticity on two-site anisotropy in the triangular lattice. (a,b) Snapshot from MC simulations of the real space spin configurations and topological charge for monolayer NiI 2 , as obtained from MC simulations with artificially modified nearest-neighbour anisotropic symmetric exchange: in closer detail, the two-site anisotropy terms have been swapped between Ni 0 -Ni 2 (5) and Ni 0 -Ni 3 (6) pairs, as shown in (c). The resulting spin-texture shown in (a) corresponds to a bi-skyrmion lattice, with vorticity m = 2 as opposed to the spontaneous anti-bi-skyrmion lattice with m = −2. Accordingly, the spin configuration of the magnetic vortices accompanying the A2Sk-topological core also changes from "half"-skyrmion-wise to "half"-antyskirmion-wise.  Fig. S3 and Fig. S4, confirm that the spin configuration (in particular its topology and chirality) is largely determined by the anisotropic symmetric exchange.

SIII. SUPPORTING MC AND DFT TESTS
FIG. S6: Skyrmionic size. As further shown in next Fig. S7, the periodicity of the skyrmionic lattice and, then, the approximate size of the a single skyrmionic object are primarly determined by the strength of the magnetic frustration, i.e. the J 3iso /J 1iso ratio. From the predicted exchange interactions in monolayer NiI 2 reported in Table I of the main manuscript, we found that the minimal size L of the supercell, that is the magnetic unit cell (m.u.c), needed to accomodate the A2Sk-kattice is 8a 0 , a 0 being the lattice constant of the unit cell. This value comes from the estimate of the propagation vector-q obtained by the evaluation of the spin structure factor S(q) (Eq. 2 in the main-text) and it is also in agreement with the minimization of the isotropic exchange interaction in momentum space J(q), given, e.g., in Ref [3]. In particular, the m.u.c. of the A2Sk lattice in NiI 2 (dashed lines inside the hexagonal 24x24 cell) comprises three anti-biskyrmions, each surrounded by six vortices. The number N of topological objects scales up with the size L of the magnetic supercell: L × L = nLm.u.c × nLm.u.c = n 2 (Lm.u.c × Lm.u.c); e.g. 24 × 24 = 3 2 (8 × 8) → N 8×8 = 3 → N 24×24 = 3 2 · 3 = 27. The radius of each anti-biskyrmion can be then estimated as r ∼ Lm.u.c/2 √ 3; the A2Sk-diameter counts 5 spins and thus it can be classified as atomic-scale skyrmionic structure [4,5].  We report here the coefficient Ai for a given nickel ion i of the SIA contribution to the magnetic interaction, which can be here simplified as HSIA = AiS 2 iz . Digits in brackets are below the accuracy of our DFT calculations (estimated to be within ∼ 10 −4 eV ). In fact, for NiI2, we also verified that no further terms appear in the expansion of the single-ion anisotropy Hamiltonian [6]: coefficients Axy, Axz, Ayz, and (Ayy − Axx) are in fact equal to zero; only the (Azz − Axx) is +0.6 meV, resulting thus into an easy-plane anisotropy. Moreover, we found that the SIA does not affect the spontaneous stabilization of the anti-biskyrmions lattice, being its contribution negligible with respect to the dominant J 1iso and J 3iso interactions and also with respect to the Jyz anisotropic term (|A/J 1iso | 0.08; |A/J 3iso | 0.10; |A/Jyz| 0.4). In fact, we performed MC calculations both without SIA and changing its sign, finding that the topology of the A2Sk and Sk lattices was unaffected.  SIII: J 1iso , J 2iso and J 3iso , and two-site anisotropy (J two−site aniso ) components (in meV) estimated for monolayer NiI2 with lattice parameter fixed to the experimental one known for the bulk phase, i.e. a = 3.89Å , within the the Liechtenstein approach (U = 1.8 eV, J = 0.8 eV ).