Interplay between Single-Ion and Two-Ion Anisotropies in Frustrated 2D Semiconductors and Tuning of Magnetic Structures Topology

The effects of competing magnetic interactions in stabilizing different spin configurations are drawing renewed attention in order to unveil emerging topological spin textures and to highlight microscopic mechanisms leading to their stabilization. The possible key role of the two-site exchange anisotropy in selecting specific helicity and vorticity of skyrmionic lattices has only recently been proposed. In this work, we explore the phase diagram of a frustrated localized magnet characterized by a two-dimensional centrosymmetric triangular lattice, focusing on the interplay between the two-ion anisotropy and the single-ion anisotropy. The effects of an external magnetic field applied perpendicularly to the magnetic layer, are also investigated. By means of Monte Carlo simulations, we find an abundance of different spin configurations, going from trivial to high-order Q skyrmionic and meronic lattices. In closer detail, we find that a dominant role is played by the two-ion over the single-ion anisotropy in determining the planar spin texture; the strength and the sign of single ion anisotropy, together with the magnitude of the magnetic field, tune the perpendicular spin components, mostly affecting the polarity (and, in turn, the topology) of the spin texture. Our analysis confirms the crucial role of the anisotropic symmetric exchange in systems with dominant short-range interactions; at the same time, we predict a rich variety of complex magnetic textures, which may arise from a fine tuning of competing anisotropic mechanisms.


Introduction
The competition between different magnetic interactions, leading to the so-called frustration, is a key ingredient for the stabilization of noncollinear and noncoplanar magnetic configurations. In particular, anisotropic interactions play a crucial role in the formation of exotic and topological spin textures, such as Bloch or Neél type skyrmions and anti-skyrmions, most commonly having a topological charge |Q| equal to one [1][2][3]. Nevertheless, topological spin textures go beyond conventional skyrmions, from high-Q (anti)skyrmionic [4][5][6][7] to halved-Q (anti)meronic spin textures and lattices [4,[8][9][10][11][12][13]. The different types of isolated two-dimensional (2D) topological spin textures are primarily characterized by the polarity p and vorticity m, whose product, Q = m · p, defines their topology. The polarity is associated with the out-of-plane magnetization profile, when moving from the core of the topological object to its edge (to infinity in a continuous description of isolated object); the vorticity is determined by the in-plane magnetization rotation and it is also referred to as the winding number, allowing only for integer values [8,14]. Therefore, (anti)skyrmions, which are characterized by reversed magnetization directions when comparing the core with its edge, corresponding to a spin configuration
Monte Carlo (MC) simulations have shown that the combination of these two properties result in the spontaneous stabilization of a topological spin structure with well defined topology and chirality.
In closer detail, the model Hamiltonian for classical spins of unit length s i = 1 is defined as: where J ij and A i are tensors describing the exchange interaction and single-ion anisotropy, respectively [53]. It is convenient to decompose the exchange coupling tensor into an isotropic part J iso ij = 1 3 TrJ ij and an anisotropic symmetric part J S ij = 1 2 (J ij + J T ij ) − J iso ij I, herein also referred to as two-ion anisotropy. Due to the inversion symmetry of the lattice (with D 3d point group), the DM interaction, which corresponds to the antisymmetric exchange term J A ij = 1 2 (J ij − J T ij ), is identically zero. Reference magnetic parameters for NiI 2 monolayer, evaluated from density-functional theory (DFT) calculations, are given in Table 1. The interactions were estimated by means of the four-state energy mapping method [53], performing PBE+U+SOC calculations [54,55] (U = 1.8 eV, J = 0.8 eV) via the VASP code [56,57]. We used a 5 × 4 × 1 supercell to estimate the SIA and first nearestneighbor interactions, and a 6 × 3 × 1 supercell for the estimate of the third nearestneighbor interaction. Supercells were built from the periodic repetition of the optimized NiI 2 monolayer unit cell, with the lattice parameter of about 3.96 Å and a vacuum distance of 17.5 Å between periodic copies of the free-standing layer along the c axis. The matrix elements of the two-ion anisotropy tensor (J S or TIA) between nearest-neighbour spins, are expressed in a local cartesian {x, y, z} basis, where x is parallel to the Ni-Ni bonding vector, and to the lattice vector a. Further details can be found in the Methods section of Ref. [40].
Taking into account such terms, the phase diagram of the magnetic system, arising when tuning the single-ion anisotropy and applied out-of-plane magnetic field, has been studied within a Monte Carlo approach. In particular, MC calculations were performed using a standard Metropolis algorithm on L × L triangular supercells with periodic boundary conditions. Starting from high temperature (T), at each simulated T, we used 10 5 MC steps for thermalization and 5 × 10 5 MC steps for statistical averaging. The lateral size of the simulation supercells was chosen as L = nL m.u.c. , where n is an integer and L m.u.c. is the lateral size of the magnetic unit cell, in units of the lattice constant a 0 , needed to accommodate the lowest-energy noncollinear helimagnetic spin configurations. Accordingly, we estimated L m.u.c. as 1/q, where q is the length of the propagation vector q minimizing the isotropic exchange interaction in the momentum space, J(q). Using the magnetic parameters listed in Table 1, the propagation vector for the isotropic model is given by q = 2 cos −1 [(1 + 1 − 2J 1iso /J 3iso )/4] [26,33], resulting in L m.u.c. 8. The results presented are thus obtained by means of calculations performed on supercells with lateral size L = 3L m.u.c. = 24.
Further insights on the magnetic configurations are obtained by evaluating the spin structure factor: where r i denotes the position of spin s i and N = L 2 is the total number of spins in the supercell used for the MC simulations. The bracket notation is used to denote the statistical average over the MC configurations. The spin structure factor provides direct information on the direction and size of the propagation vectors. The topological nature of the multipleq phases has been assessed by evaluating the topological charge (or skyrmion number) of the lattice spin field of each supercell as Q = ∑ i Ω i ; Ω i is calculated for each triangular plaquette as [58] tan 1 2 In the following, starting from the exchange interactions estimated in the prototypical monolayer NiI 2 as representative example of highly frustrated semiconducting 2D triangular lattice systems, we explore the magnetic phase diagram as a function of the strength and direction of the SIA as well as of an applied out-of-plane magnetic field B z .

Results
In Figure 1, we show the low temperature (T = 1 K) magnetic phase diagram in the SIA-field (A zz −B z ) plane, obtained by fixing the exchange coupling interactions to the reference values reported in Table 1, and varying the value and direction of the single-ion anisotropy; for each A zz value, an increasing magnetic field (B z ) was also applied along the perpendicular direction of the monolayer triangular lattice. As shown in Figure 1a, for each combination of exchange coupling and A zz parameters, there are field-induced topological phase transitions between states with different topological charge Q, as listed below, approximately, per magnetic unit cell (m.u.c.): As schematically represented in Figure 1c, the phase diagram consists of various trivial and topologically equivalent phases-i.e., displaying the same Q m.u.c -but characterized by different kinds of spin configurations. As shown in Figure 1b, for different realizations of the easy-axis and easy-plane anisotropy, field-induced phase transitions are accompanied by discontinuities of the out-of-plane magnetization M z . Within the considered range of applied B z , the magnetization saturation M z = 1, corresponding to all spins aligned parallel to the magnetic field, is reached only in the case of easy-axis anisotropy and for small easy-plane anisotropy; indeed, the latter penalizes the out-of-plane component of the spins. The in-plane magnetization is zero in all cases.
In the following, we describe and analyze the various anisotropy-field induced magnetic orders (Figures 2-6), also providing visual details of the out-of-plane and in-plane components of spins for the various phases (Figures 7 and 8  Evolution of the out-of-plane magnetization M z as a function of the magnetic field B z for increasing SIA, easy-axis (A zz /|J| 1iso < 0, top) and easy-plane (A zz /|J| 1iso > 0, bottom), respectively. (c) Schematic A zz − B z phase diagram identifying regions occupied by the different spin configurations, labelled with capital letters. Phase boundaries must be regarded only as semi-quantitative, as the field-induced topological phase transitions go through almost continuous transformations of the spin texture and metastable spin states with fractionalized topological charge, preventing the identification of accurate and sharp phase boundaries. On the right-side of the y-axis, A zz /|J yz | ratio is also reported. The horizontal solid line separates the easy-axis (above) and easy-plane (below) SIA regions.

Topologically Trivial Spin Orderings
We start our analysis by discussing the topologically trivial spin configurations found in the strong-field region of the phase diagram, and in the weak-field region for strong SIA (both easy-axis and easy-plane). In Figure 2, we show both spin configurations and topological charge densities (Ω i ) alongside the associated spin structure factor of such phases, labeled as (A), (B) and (H) in the phase diagram of Figure 1c. Figure 2. Trivial spin states with Q m.u.c = 0. Snapshots at T = 1 K from Monte-Carlo simulations of real-space spin configurations (a,d,g,j) and corresponding topological charge densities Ω i (b,e,h,k), in sequence: black arrows represent in-plane {s x , s y } components of spins; colormap indicates the out-of-plane s z spins component in the first snapshot and Ω i in the second one. The associated spin structure factor S(q) is also shown (c,f,i,l); its decomposition in the Cartesian components S x (q), S y (q) and S z (q) is shown for the spin spirals configurations (first and second panels). Capital letters on each spin configuration refer to the labels used in the phase diagram Figure 1c to name each spin configuration: A-phase refers to a single-(q) spiral (proper screw), B-phase refers to a single-(q) spiral (tilted cycloid), H-phase refers to a triple-(q) trivial state with FM component at Γ-(q = 0), I-phase refers to a FM state.
For weak fields, two helical single-q states (A and B phases) appear for strong easyaxis and easy-plane SIA, whose spin configurations are shown in Figure 2a,d, respectively. No finite topological charge appears in such phases, as clearly shown in Figure 2b,e. From Figure 2c,f, the propagation vector is q 1 = (−2δ, δ), with q/2 = δ 1/8 in reduced coordinates, being mostly determined by frustrated isotropic exchange. Due to the threefold rotational symmetry of the triangular lattice, the single-q helices propagating along the x-axis with q 1 are equivalent (and, hence, energetically degenerate) to helices propagating along the symmetry-equivalent directions rotated by +120 • , with q 2 = (δ, −2δ), and −120 • , with q 3 = (δ, δ). The plane of spins rotation is selected by both SIA and TIA, which thus determine the nature of the helical states. The inspection of both spin structures and spin structure factors reveals that: (i) phase (A) consists of a proper-screw spiral propagating along the Cartesian x-axis, i.e., along the Ni-Ni bond direction, with spins rotating in the perpendicular yz-plane; (ii) phase (B) is a tilted cycloid, i.e., a helix where spins rotate in a plane containing the propagation vector (parallel to x) but tilted around it, causing the spins to acquire a non-zero z-component highlighted by the peaks of S z (q) (cfr Figure 2f). The modeling of the single-q proper-screw and cycloidal helices is presented in Appendix A.
The combination of the three spirals defines a triple-q state, that is the case for the H-phase depicted in Figure 2g. This is a trivial noncollinear spin configuration, but characterized by local nonzero scalar chirality, as evident from the colormap of the topological charge density in Figure 2. Indeed, it is possible to recognize a hexagonal lattice formed by six vortices (m = +1) with chirality opposite to the surrounded anti-bi-vortex core (m = −2). This state occupies a wide region of the phase diagram for high values of applied magnetic field, before evolving into a pure ferromagnetic state [phase (I), Figure 2j]. In addition to the weak peaks at q 1 , q 2 , and q 3 (Figure 2i), this H-phase displays an out-of-plane ferromagnetic component, as shown by the relevant peak in the S(q) at the Γ-point (q = 0), which is also reflected in the remarkable out-of-plane component of the magnetization, 0.5 M z < 1 (Figure 1b, as well as in the s z profile depicted in Figure 8m). Noteworthy, such triple-q trivial state is largely favored by the easy-plane anisotropy, which indeed competes with B z , penalizing the out-of-plane spin component; indeed, for A zz /|J| 1iso > 0.15 (A zz /|J yz | > 0.75), the FM state is never achieved within the wide range of explored B z . Conversely, the easy-axis anisotropy sustains the B z action, favoring the out-of-plane spin orientation; accordingly, the region of the phase diagram occupied by the H-phase reduces, until it disappears for A zz /|J| 1iso −0.75 (Figure 1c).

Topological Spin Orderings
In this section, we turn our attention to the various topological spin configurations which can be realized by tuning the competition between the single-ion anisotropy and the applied field. The related spin textures, topological charge densities and spin structure factors are shown in Figures 3-5, respectively. All topological lattices are triple-q states, Figures 3c,f, 4c,f and 5c], characterized by atomic scale skyrmionic or meronic structures composed by nano-sized topological objects with a radius counting 4 spins, and thus a diameter of few units of the lattice parameter a 0 .
In Figure 3a,d, we show the high-Q topological lattices labelled C-and D-phases, respectively, which occupy the region of the phase diagram for zero and low magnetic field; in particular, the C-phase is the spontaneous topological phase reported in [40], and stabilized by the exchange magnetic couplings given in Table 1. In closer detail, the topological lattice of the C-phase ranges over −0.30 A zz /|J| 1iso +0.50 and −1.50 A zz /|J yz | +2.50 for zero field, whereas it slightly extends over the easy-plane anisotropy part of the phase diagram upon low B z , as shown in Figure 1a,c. It is a triple-q state characterized by a hexagonal lattice formed by six vortices (V, m = +1) with downward central spin, hosting at the center of each hexagon an anti-bi-vortex (A2V, m = −2) with upward core, as shown in Figure 3a. Such spin texture defines a homogeneous topological charge density, as depicted in Figure 3b, hence uniform scalar chirality, and gives rise to a topological charge of six per magnetic unit cell. This topological lattice can be regarded as the periodic repetition of topological objects, consisting of two vortices and one anti-bi-vortex (schematically depicted in Figure 7g), all contributing to defining a global |Q| = 2. As evident from the s z profile reported in Figure 7a, the central A2V does not display a uniform, unitary polarity along its surrounding perimeter (highlighted as a purple circle traced in Figure 7b): the upward spin of the core (s z = +1) is not reversed downward (s z = −1) at all points of its finite edge. Therefore, it brings a fractionalized topological charge Q A2V = −2 · p A2V , with 0 < p A2V < +1. The missing Q fraction is carried by the two vortices, which, similarly, bring a fractionalized charge Q V = +1 · p V , with −1 < p V < 0. Accordingly, |Q| = |Q A2V + 2Q V | = 2. The minimal magnetic cell accommodating this topological lattice consists of the composition of three repeated objects; hence |Q| m.u.c = 6. The topological lattice can therefore be interpreted as a fractionalized anti-bi-skyrmion (A2SK) lattice, where the fractionalization of the topological charge can be ascribed to the A2SKs close packing, leading to the incomplete spin wrapping highlighted in Figure 7 [10]. A similar realization of a fractionalized skyrmion lattice has been recently reported in MnSc 2 S 4 , where an applied field has been experimentally shown to stabilize a lattice of fractionalized Bloch-type skyrmions and incipient merons [59].
By increasing the easy-axis anisotropy, i.e., for −0.75 A zz /|J| 1iso −0.45 and −3.75 A zz /|J yz | −2.25, keeping low values of the applied magnetic field, the magnetic phase transforms into the D-phase shown in Figure 3d. Despite the weak FM component seen in the spin structure factor (Figure 3f), which is compatible with both applied field and enhanced easy-axis anisotropy, the spin configuration approaches an ideal triangular lattice of anti-bi-skyrmions: each anti-bi-vortex is surrounded by a magnetic background, where almost all spins are fully reversed with respect to the core of the A2V along the radial directions, as depicted in Figure 8a,b. Nevertheless, the proximity of the small sized A2SKs still produces a weak overlap between these topological objects, in turn, causing the formation of residual vortices at the center of the triangles formed by three nearest-neighbor A2SKs [10]; therefore, these anti-bi-skyrmions carry a fractionalized topological charge, that is 1.5 |Q A2SK | < 2 (or equivalently 4.5 |Q| m.u.c < 6 considering three anti-biskyrmions per magnetic unit cell). Accordingly, the topological charge density map is no longer homogeneous, rather the highest intensities are localized around the anti-bi-vortex core (Figure 3e). We finally notice that the cores of each A2SK is found to align anti-parallel to the applied field direction. The spin texture shown in Figure 3d has been obtained by applying a negative B z to ease the comparison with the C-phase displayed in Figure 3a; we verified that a reversal of the magnetic field direction systematically causes a reversal of the A2SK-core magnetization and an alignment of the magnetic background with the field.

Topological
Lattice with |Q| m.u.c = 3 In Figure 4a,d, we show the topologically equivalent |Q| m.u.c = 3 lattices, labelled E-and F-phases, respectively, which occupy a wide portion of the phase diagram in an intermediate range of applied magnetic field B z (Figure 1c). Despite carrying the same topological charge, their spin textures are markedly different.  As shown in Figure 4a, the E-phase consists of a triangular skyrmion lattice, with vortices and anti-bi-vortices appearing at the intersection of clearly visible Bloch-like skyrmions, namely, at the center of the triangles that build up the skyrmion lattice. Such phase can be seen as a direct transition of the stable C-phase induced by intermediate values of the applied perpendicular field: upon B z , spins of the anti-bi-vortex core remain parallel to the field, spins of one of the two vortices also tend to align with the magnetic field, while spins at the core of the other vortex remain anti-parallel. Such a change in the s z spins component in one of the two vortices affects the spin scalar chirality, thus leading to the topological transition from |Q| m.u.c = 6 in the C-phase to |Q| m.u.c = 3 in the E-phase. In the latter phase, the topological charge density map shows in fact negligible contributions from the A2V, while giving rise to an opposite sign of Ω i for the upward and downward vortices. The highest contribution to the topological charge arises from the vortices whose core-magnetization remains antiparallel to the applied field, thus defining the Bloch-type skyrmion lattice. Indeed, the corresponding s z profile, shown in Figure 8g,h, testifies for an almost full reversal of the out-of-plane spin component, i.e., polar angle variation of about 180 • , when moving from the core (s z = −1) to the edge (s z = +1) of the Bloch-type vortex.
Such skyrmion lattice is found to be robust with respect to easy-axis anisotropy. On the other hand, a strong enough single-ion easy-plane anisotropy (A zz >> 0) induces its transformation into the F-phase, observed in a narrow region of the investigated phase diagram (see Figure 1c): the out-of-plane component of magnetization antiparallel to the applied field is here further suppressed by the strong easy-plane SIA. In particular, the analysis of the spin texture shown in Figure 4d, combined with the s z profile (Figure 8d,e) and the topological charge, indicates that the F-phase is a lattice of anti-bi-merons (A2M): the A2V of the C-phase preserves the upward core, with spins parallel to B z , while the surrounding vortices preserve only the planar spins components. Indeed, the polar angle changes by 90 • when moving along the radial direction from the core to the edge, i.e., p = +0.5, as from the s z profile in Figure 8d. As a consequence, the topological charge of the F-phase is halved with respect to that of the C-phase, that is |Q| m.u.c = 6 → 3. To the best of our knowledge, this is the first time that such kind of A2M texture is reported.

Topological
Lattice with |Q| m.u.c = 1.5 As a last relevant topological phase composing the phase diagram, we show in Figure 5a a meronic lattice with |Q| m.u.c = 1.5, which can be regarded as an intermediate spin configuration during the evolution of the SK-lattice of the E-phase into the trivial triple-(q) ferromagnetic H-phase upon B z (Figure 1c). The competition of the intermediate-strong easy-plane anisotropy with the intermediate-strong applied field causes the transformation of the Bloch-like skyrmion into a meron with halved topological charge. As shown by the spin texture in Figure 5a and the s z profile of the vortex core in Figure 8j,k, the latter loses the spin component along the z-direction, while it remains aligned to the field in the surrounding vortices and anti-bi-vortices, which occupy alternatively the center of the triangles formed by the merons: spins at the edge of the topological vortex are directed perpendicularly with respect to its core, halving the topological charge of the Bloch-like skyrmion lattice of the E-phase, |Q| m.u.c = 3 → 1.5.  Figure 3, snapshots at T = 1 K of the real-space spin configuration (a) and topological charge density Ω i (b), for the G-phase, i.e., a triple-(q) state, consisting of a lattice of merons from the action of large easy-plane anisotropy and B z , as from the q = 0-peak in the S(q) (c), on the SK-lattice of the E-phase.

Field-Induced Phase Transition from Two-Ion Anisotropy Tuning
In Figure 6, we report the evolution of the topological charge per magnetic unit cell |Q| m.u.c as a function of increasing B z , and for different values of the off-diagonal J yz exchange coupling term, which is treated here as a measure of the non-coplanarity induced by the two-ion anisotropy. SIA is kept fixed to the NiI 2 reference value (Table 1); in fact, from the phase diagram previously discussed (Figure 1), we found that the SIA is basically ineffective against the anisotropic symmetric exchange interaction for a wide range of values (i.e., −1.50 A zz /|J yz | +2.50 and −0.30 A zz /|J| 1iso +0.50).  Figure 6. Evolution of |Q| per magnetic unit cell (L m.u.c. = 8a 0 ), as a function of the magnetic field B z for different off-diagonal terms, while keeping fixed diagonal exchange parameters and easy-plane anisotropy to values reported in Table 1. The color gradient of the solid points evolves with the A zz /|J yz | ratio. At zero and low B z field, the system stabilizes a proper screw like the A-phase for J yz /J 1iso ≤ 0.14, and the topological lattice of the C-phase for J yz /J 1iso > 0.14. At intermediate B z , the system undergoes a topological phase transition to a Bloch-type skyrmion lattice, E-phase. At high B z , the transition to the ferromagnetic state, I-phase, takes place passing trough the trivial triple-q state, H-phase.
For J yz > 0.14 J 1iso (the reference value is J yz 0.20 J 1iso ), the topological lattice with |Q| m.u.c = 6, represented by the C-phase (Figure 3a), is the lowest energy spin configuration stabilized by the considered magnetic interactions. In particular, in line with the results shown in Figure 1a, two sharp topological phase-transitions, signalled by the abrupt change of the total topological charge, are induced under applied magnetic fields: |Q| m.u.c = 6 → 3 → 0. The |Q| m.u.c = 3-state obtained under intermediate B z is the Bloch-type skyrmion lattice of the E-phase; for strong B z , the transition to the trivial ferromagnetic I-phase takes place passing through the intermediate H-phase, characterized by the topologically trivial triple-q state (Figure 2j,g, respectively).
For J yz ≤ 0.14J 1iso , two topological phase transitions still take place, but the stable phase at zero and weak magnetic field is now a trivial single-q state; |Q m.u.c | = 0 → 3 → 0 as a function of B z . In particular, such single-q state is found to be a proper-screw spinspiral. The magnetic phases induced by the applied fields are closely related to those observed when J yz > 0.14J 1iso . Therefore, a Bloch skyrmion lattice can still be stabilized by a magnetic field applied on a topologically trivial helical state, when the non-collinearity brought in by the two-ion anisotropy is not sufficiently strong to drive a spontaneous high-Q topological lattice.

Discussion
In the previous section, we have described the rich anisotropy-field phase diagram (Figure 1) obtained by tuning strength and direction of the single-ion anisotropy (from easy-plane to easy-axis) as well as upon increasing magnetic field, while keeping fixed the strength of both the isotropic and anisotropic symmetric exchange interactions describing the short-range magnetic interactions of the spin-spin Hamiltonian (1) for a centrosymmetric triangular lattice. We found that, considering the magnetic interactions reported in Table 1, the spontaneous C-phase, consisting of a hexagonal lattice of vortices (m = +1) and anti-bi-vortices (m = −2) with a total topological charge of six per magnetic unit cell, is a persistent thermodynamically-stable phase in a wide portion of the SIA-field phase diagram; it is, to a wide extent, independent on the easy-plane or easy-axis character of the single-ion anisotropy. Our analysis also unveils different field-induced topological phase transitions that can be realized for different values of the SIA coupling constant.
Nevertheless, it is important to notice that all the different emerging triple-(q) states, ranging from topological (phases C, D, E, F, G) to trivial (H-phase) spin configurations, exhibit a common planar spin texture. To better appreciate this feature, we show in Figures 7c and 8c,f,i,l,o the color-gradient plots of the in-plane spin components highlighting the azimuthal angle θ ∈ [0 • , 360 • ), defined by the s x and s y components of the spin vector s, and which is related to the in-plane spins rotational direction.  (Table 1). Spins on the nearest-neighbour M sites of the central site orient according to the in-plane projection of the noncoplanar principal axes, as from the zoom on the spin texture of the A2V-core obtained by the MC simulations (f).
Such specific in-plane orientation of spins can be traced back to the TIA, which tends to orient the spins along given noncollinear and noncoplanar orientations in space, defining both the helicity and vorticity of the spin pattern and thus behaving as an emergent chiral interaction. Indeed, as discussed in Ref. [40], the principal magnetic axes per bond of the J S tensor are not all parallel either to the lattice vectors or with the principal axes of neighbouring bonds, as shown in Figure 7d. This introduces noncoplanar components in the spin-spin interaction and, thus, frustration in the relative orientation of the spins, which ultimately arises from the non-coplanarity of the spin-ligand M − X plaquettes mediating the exchange interactions in the spin triangular lattice. In particular, the in-plane projection of the principal axis directed along the ligand X − X direction [indicated by red arrows in Figure 7d,e] fixes the in-plane orientation of the nearest-neighbour magnetic moments giving rise to the A2V spin-pattern [ Figure 7e,f]. The accommodation of the anti-bi-vortices in the spin lattice, which tend to overlap because of the short-period modulation due to the isotropic magnetic frustration (J 3iso vs J 1iso ), causes the emergence of the surrounding vortices. With the given strong anisotropic symmetric exchange interaction considered in the present case, the two-ion anisotropy dominates over the single-ion anisotropy and external magnetic field: SIA and B z tune the out-of-plane spins component, modifying the polarity p, and thus the final spins scalar chirality and related topology (Q) of the spin configuration; the planar spin texture results always fixed by the TIA. All the triple-(q) phases can be in fact seen as a transformation of the initial C-phase via major modifications of the s z component of spins; the in-plane directions of the spins remaining unchanged. Moreover, even though single-q states are found to be energetically more stable than triple-q ones when A zz 2.5|J yz | and A zz −1.5|J yz | (or also when J yz ≤ 0.14J 1iso ), an applied out-ofplane magnetic field can still induce a transition to a topological skyrmion lattice, whose in-plane spin configuration appears to be determined by the two-ion anisotropy.
Beyond the fundamental interest of our findings, which identifies the distinct role played by different competing magnetic interactions, the rich anisotropy-field phase diagram calls for further studies. In particular, the investigation of magnetic properties in other two-dimensional van der Waals materials and their possible interplay with dielectric properties, as also pointed out by recent experimental achievements [60,61], open the path to new low-dimensional magnetoelectrics and/or multiferroics.

Conclusions
In this work, we have theoretically investigated the effects of competing single-ion and two-ion anisotropies in a triangular lattice with strong magnetic frustration, as occurring in monolayers of van der Waals nickel dihalides. By means of Monte Carlo calculations, we analyzed the parameter space spanned by SIA and applied B z field for a spin-lattice model whose parameters have been estimated for a monolayer of prototypical semiconducting NiI 2 . Our analysis has revealed a rich phase diagram comprising different magnetic phases, from topologically trivial single-(q) and triple-(q) states to topological triple-(q) states. The strong magnetic frustration, arising from the competing isotropic nearest-neighbour FM and third nearest-neighbour AFM exchange interactions, promotes the onset of shortperiod helimagnetic configurations, whereas the strong two-ion anisotropy combined with the geometrical frustration of the underlying triangular lattice favours the stabilization of triple-q states. At zero magnetic field, these result in a lattice of vortices and anti-(bi)vortices carrying a total topological charge per magnetic unit cell |Q| m.u.c = 6, which is robust within a wide range of single-ion anisotropy strength. Such topological lattice can be interpreted as a fractionalized anti-bi-skyrmion lattice where each anti-bi-skyrmion is surrounded by six vortices arising from the overlap with neighbouring anti-bi-skyrmions, causing a fractionalization of the topological charge of individual anti-bi-vortices. Both the single-ion anisotropy and the applied field act primarily on the out-of-plane component of the spin configurations, hence on the polarity p of the vortical states, modulating both the size of anti-bi-vortices and the fractionalization of their charges. As a general trend, on the one hand, easy-axis anisotropy is found to increase the localization of anti-bi-vortices, thus reducing the "spilling" of topological charge to surrounding vortices; on the other hand, a strong easy-plane anisotropy, when combined with applied field, would eventually favour a crossover to meronic lattices, i.e., lattices formed by topological objects with half-integer polarity. Nevertheless, an out-of-plane magnetic field, sustained by the strong anisotropic exchange, is always found to trigger a topological transition to a Bloch-like skyrmion lattice, when applied either on a single-or a triple-q state. Interestingly, the inplane components of the spin texture are found to be extremely robust across almost the whole phase space explored in this work, with exceptions only for extreme values of SIA and fields.
In conclusion, the different topological phases forming the SIA-field phase diagram can be understood by a change of their polarity, which is directly tuned by both the singleion anisotropy and the applied field, whereas their vorticity can be always traced back to the two-ion anisotropy on the triangular lattice, which is ultimately responsible for the promotion of both topological and trivial triple-q states. Therefore, our study corroborates the potential role of frustrated anisotropic symmetric exchange in defining the vorticity of noncollinear and noncoplanar spin configurations, leading to various possible topological lattices in 2D semiconducting magnets.
Author Contributions: D.A. performed Monte Carlo simulations and analyzed the results. P.B. developed the code for Monte Carlo simulations, for the calculation of the topological charge Q and for the spin structure factor S(q), and contributed to the analysis. D.A. and P.B. wrote the manuscript. S.P. supervised the project. All the authors discussed the results. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: Main results are reported in this article. All other data that support the findings discussed in this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The single-q proper-screw (ps)-phase (A)-and cycloidal (c)-phase (B)-helices propagating along the x axis can be modeled, respectively, as where α i = qx i and θ is the polar angle describing the tilting of the spin-rotation plane in the cycloidal configuration. The corresponding mean-field energies evaluated for model Equation (1) can be generally written as: where E iso (q) = 2J 1iso cos q 2 + (J 1iso + 2J 3iso ) cos q + J 3iso cos 2q, for all single-q states, while E SI A ps = A zz /2 for the proper-screw spiral and E SI A c = A zz cos 2 θ/2 for the cycloid. The mean-field expression for TIA contribution is less transparent, but one can easily show that the off-diagonal TIA term J S yz only affects the cycloid meanfield energy contributing as J S yz sin 2θ (cos q/2 − cos q). The interplay of both SIA and TIA removes the degeneracy of all coplanar helical states arising from isotropic magnetic frustration. A strong easy-axis SIA (A zz < 0) would promote a proper-screw or a vertical (θ = 0) cycloidal helix, with spins rotating in the yz or xz plane, respectively. The degeneracy between these two helices is then removed by TIA, being E TI A ps − E TI A c = J ⊥ (cos q/2 + cos q), where 2J ⊥ = J S xx − J S yy . Since J ⊥ < 0 from Table 1, the proper-screw helix of phase (A) is ultimately stabilized over the vertical cycloid by the two-ion anisotropy. In the opposite limit of strong easy-plane SIA (A zz > 0), the in-plane θ = π/2 cycloid is energetically more favourable than the proper-screw spiral, while the tilting of the spin-rotation plane with respect to the xy plane is induced by the strong off-diagonal J S yz , pointing to TIA as the microscopic origin of the B-phase.