Theoretical investigation of antiferromagnetic skyrmions in a triangular monolayer

The chiral spin textures of a two-dimensional (2D) triangular system, where both antiferromagnetic (AF) Heisenberg exchange and chiral Dzyaloshinsky-Moriya interactions co-exist, are investigated numerically with an optimized quantum Monte Carlo method based on mean-field theory. We find that: helical, skyrmionic and vortical AF crystals can be formed when an external magnetic field is applied perpendicular to the 2D monolayer; the sizes of these skyrmions and vortices change abruptly at several critical points of the external magnetic field; each of these AF crystals can be decomposed into three periodical ferromagnetic (FM) sublattices. The quantum ingredient implemented into the theoeretical framework helps to track the existence of AF skyrmion lattices down to low temperatures.


I. INTRODUCTION
Skyrmions in ferromagnetic materials have been intensely studied in recent years. [1][2][3][4][5][6][7][8][9][10] Their small size, in a range of 1 nm ∼ 100 nm, and the very weak electronic current, around 10 6 Am −2 , required to drive them to motion, 11 make them ideal candidates in future data storage and racetrack memory devices. Unfortunately, the Magnus force, 12-14 that acts transversely on the ferromagnetic (FM) skyrmions by the applied electric current, eventually pushes them off the nano-track edge, which greatly hinders their applications in the aforementioned memory devices. 15 In contrast, an antiferromagnetic (AF) skyrmion in a square lattice can be decomposed into two almost identical FM sublattices. Thus, the Magnus forces acting on the two FM sublattices can be completely cancelled, 16 so that the AF skyrmion can move much faster straightly along the direction of the applied electric current. 16,17 AF skyrmions have been investigated theoretically in 2D square and triangular antiferromagnets. [18][19][20] In their Monte Carlo study, Keesman and his colleagues discovered that an isolated skyrmion can be stabilized by a strong external magnetic field at nonzero temperatures, but only in a very tiny finite-sized 8 × 8 square antiferromagnet. 20 For triangular lattices, Rosales et al. 19 observed with classical Monte Carlo (CMC) simulations that when AF Heisenberg exchange (HE) and Dzyaloshinsky-Moriya (DM) interactions co-exist, AF skyrmionic lattices can be induced by an external magnetic field applied normal to the lattice plane in a very low temperature range. However, they found no evidence of AF skymrion lattice (SL) formation in the AF chiral square lattice. The latter aspect was addressed by one of us 28 utilizing a quantum simulative method, namely the SCA approach, where a self-consistent algorithm is used in the frame of quantum theory. [21][22][23][24][25][26][27][28][29] It was found that both Bloch and Néel types AF SLs could be induced in a 2D lattice of square crystal structure by a very strong external magnetic field exerted normal to the monolayer plane. 28 We expect that quantum effects, missing in CMC simulations, to be important in AF materials, which motivates the current study by reinvestigating AF SLs formed in the 2D triangular lattices with a method built upon quantum theory.
We utilize here an optimized quantum Monte Carlo (OQMC) method based on a meanfield approach. The Metropolis algorithm is employed, and a simple trick is adopted to improve the computational efficiency. Astonishingly, by using this method we are able to plot, for examples, the well periodic and symmetric AF vortical lattices (VLs), AF skyrmionic lattices (SLs) of both Bloch and Néel types, which can appear at any temperatures in 2D AF chiral magnets, just with the computational results obtained from the last iteration. In contrast, if CMC method is applied to calculate a spin configuration at a finite temperature, averages must be made over thousands of last iterations.
Using the OQMC method, we simulate and investigate in detail the chiral spin textures of a 2D triangular system where the AF HE and chiral DM interactions are both present.
Consequently, we find that, within an external magnetic field applied perpendicular to the 2D monolayer, AF HLs (helical lattices), SLs and VLs can be induced at much higher temperatures than predicted by previous authors with CMC simulations; 19 the SL states prevail in a broad T − H phase area; the sizes of the skyrmions and vortices formed in the spin crystals or lattices changes abruptly when modifying the external magnetic field; each of these AF SL and VL can be decomposed into three FM sublattices that are in fact FM SLs or FM VL, respectively. To check the accuracy of the results obtained with the OQMC method, we carry out simulations using the SCA approach for a particular case with a fixed field strength, and find the resulting AF skyrmionic lattices from both theoretical approaches to be identical.

II. THE QUANTUM MODEL AND COMPUTATIONAL ALGORITHM
The two-dimensional triangular spin-lattice is considered to be in the xy-plane, and its corresponding Hamiltonian can be expressed as Here, the first two terms represent the HE and DM interactions between a pair of spins at the i-th and j−th lattice sites, respectively. In practice, we restrict the interactions to the nearest neighboring spins. The third term stands for the uniaxial magnetic anisotropy assumed to be normal to the 2D monolayer, which is usually termed as the perpendicular magnetic anisotropy (PMA), and the last one denotes the Zeeman energy of the magnetic system placed in an external magnetic field. The strengths of the two-spin HE, DM and single-spin PMA interactions are J ij , D ij , K A , respectively. If D ij is along the r ij direction (the vector r ij connecting both spin sites), the induced skyrmion is of Bloch-type; whereas if D ij lies in the film plane and is perpendicular to r ij , the skyrmion is of Néel-type.
In light of quantum theory, the spins appearing in Eq.(1) are quantum operators, 21-29 rather than classical vectors. When S = 1 as it is assumed in the present work, the matrices of the three spin components are: respectively, in the Heisenberg representation. The thermal expectation value of a physical observable A at temperature T can be evaluated with whereÂ is the operator of the observable A, and β = −1/k B T .
When the quantum Monte Carlo (QMC) method was proposed, Metropolis algorithm was also conventionally employed. 23,27 That is, in every simulation step, a spin S i is randomly selected from the considered magnetic system, then rotated randomly within a narrow spatial cone. Afterwards, a random number r is generated to compare with p = exp(−∆E i /k B T ), where ∆E i is the energy change caused by the rotation. If r ≤ p, the operation is accepted, otherwise discarded. All simulations are started from random magnetic configurations above the magnetic transition temperatures, then carried out stepwise down to very low temperatures with a reducing step ∆T < 0.
A simple trick has been taken to optimize the QMC method: Once the rotated state of S i is accepted, the states of its neighboring spins have to be updated immediately, since their states are also affected by the operation. Otherwise, errors in the total energy of the magnetic system will be accumulated, so that the computational program may converge to an incorrect spin configuration.
The above technique taken to optimize QMC method seems very easy and simple. However, it has been proved to be very effective, leading to quick convergence of the simulations. More specifically, OMQC method makes it possible for the computational results obtained in the final iteration after convergency to accurately describe the complicated and detailed spin textures of considered magnetic system at all temperatures. For instance, with such results we are able to depict well symmetric and periodic FM and AF SkLs of the Bloch-and Néel-types in a broad T − H phase area in 2D magnetic systems.

III. COMPUTATIONAL RESULTS
Each lattice site of the AF monolayer is considered to be occupied by an S = 1 spin, both HE and DM interactions are limited to the nearest neighboring spins, and their strengths assigned to

A. Phase Diagram
From our simulated results, we find that well symmetric and periodical AF SLs can be generated when the external magnetic field falls in a wide range 1 ≤ H ≤ 7.4, as shown in the phase diagram plotted in Figure 1(a). To get this diagram, the applied magnetic field strength is varied from 0.4 to 7.8; and in an external field of fixed strength, the simulation is initiated from a random spin configuration in the paramagnetic phase, then temperature is lowered with a reducing step ∆T = -0.1 or -0.05.
In Figure 1(a), T HL , T SL , T SV represent the transition temperatures when the magnetic system condenses to the AF helical, skyrmionic and vortical lattices, respectively. Therefore, the area below the curve T α and above the curve T β is in the α phase; whereas the region below the lowest curve T γ is in the γ phase. For example, as 6.0 ≤ H ≤ 7.4, the area between the T V L and T SL curves is the VL phase, whereas the region below the T SL curve is in the SL phase; while in the field range 1.0 ≤ H ≤ 1.6, the area between the T SL and T HL curves is the SL phase, whereas the region below the T HL curve is the HL phase. In a wide field range 1.8 ≤ H ≤ 5.8, only AF SLs appear below the T SL curve; while for 0.4 ≤ H ≤ 0.8, only HL states prevail below the T HL curve.
As expected, at high temperatures the system is completely polarized by the applied magnetic field and becomes ferromagnetic. Above T = 2.8, the spins are completely polarized to the z-direction by the external magnetic field. When temperature drops to T = 2.6 and 2.5, the positive and negative branches of S x and S y , though very weak, are comparable in magnitude, thus AF SLs can be generated. However, as temperature drops further, the magnitudes of S x and S y differ considerably, while those of S y and S z are comparable, so that vertical AF HL state dominates the whole low temperature 6 range.
In the second case shown in Figure 1(c), H =7. As T > 1.6, only S z is sizable. That is, the whole system is almost completely polarized by this strong external magnetic field. Below T = 1.6, the positive and negative branches of both S x and S y grow gradually with decreasing temperature, and later become comparable in magnitude, but they are shifted relatively with each other in the vertical direction, so that VLs can be formed within temperature rangle 1.55 ≥ T > 1.3, whereas SLs are induced below T SL = 1.3.
It is easy to understand that, when H = 7.6 as shown in Figure 1(d), only VL can be observed below T V L =1.15, since all spins have been rotated by the external magnetic field toward the zdirection. Especially, when H is increased to 7.6, the size of each vortex suddenly grows, and we have to use a 60 × 60 lattice to generate these VL structures.

B. FM Sublattices of Helical, Skyrmionic and Vortical Lattices
The co-presence of AF HE and DM interactions gives rise to very complicated spin configurations, so that the AF spin lattices, though they are well periodical and symmetric, look quite puzzling if the xy projection and z-contour of such a spin crystal are plotted together. Fortunately, each of these spin textures can be decomposed into three FM sublattices.  textures in these interstitial areas seem quite 'random' in the xy-plane, but they are astonishingly still periodic in each FM sublattice. These 'disordered' textures look completely different from those obtained in CMC simulations, 19 and as a result, give rise to increased magnetic entropy, S M , and reduced total free energy, F , of the whole system, since F = E − T S M , so that the AF SLs are better stabilized.
When H = 7.6 as shown in Figure 4, the z-components of all spins become positive, so only VL can be formed below T V L = 1.15. The AF vortices form regular hexagonal VL, whose three FM sublattices show the same pattern but all rotated around the z-axis. More interestingly, the size of each vortex increases suddenly, so that we have to use a 60 × 60 lattice to produce the VL texture.

C. Distribution of Topological Charge Density of Each Sublattice
For a discretized model, the topological charge can be expressed with 19,30,31 where A L, r i = S r i · S ra × S r b denotes the local chirality, and r i , r 1 ∼ r 4 are the sites involved in the calculation of χ Q as described in Rosales et al.'s article. 19 Using above formulas, we find that the calculated topological charge density of each sublattice also forms periodical crystals, as shown in Figures 6 ∼ 7(b,c,d,), that are obtained at H = 7, T = 0.05, and H = 5, T = 2, respectively. Here, the rainbow-like contour denotes the topological charge density, ρ, at the lattice sites, and the arrows show the xy projections of the FM spin sublattice of the AF SL. Clearly, each charge density crystal is almost exactly identical with the corresponding FM spin sublattice. The three ρ sublattices of an AF SL differ from each other, and we observe much densely packed ρ unit cells in Figure 7(b,c,d).
In Figure 6(a) and Figure 7

D. Spin Structure Factors of the Helical, Skyrmionic and Vortical Lattices
To characterize the spin-textures, the static spin structure factors in the reciprocal space have to be calculated, of which the out and in plane components S ⊥ (q) and S // (q) are defined as respectively.

IV. COMPARISON WITH SCA METHOD
From now on, other simulated results will be displayed in the appendixes for compactness.
To check our results, we repeat the simulations by means of the SCA method. [21][22][23][24][26][27][28][29] Astonishingly, as H = 7, the xy-projections of the AF SLs obtained at different temperatures with the two methods, when plotted together, overlap with each other exactly, except for T = T V L = 1.55, where the system starts to condense to the VL state, only at a few sites the spins align in slightly different directions. Figure S1 shows an example at a very low temperature, T = 0.05, in comparison with the results, that are shown in Figure 3 and calculated with the same set of parameters using the OQMC method. No doubt, if the A, B and C FM sublattices obtained with the two quantum approaches are plot together correspondingly, every pair will overlap with each other almost exactly.
For the same purpose, the calculated total topological charges versus varying temperature as H = 7, and charge density textures of the A, B and C FM sublattices of the AF SL obtained at T = 0.05 with the SCA approach, are displayed in Figure S2. When compared with the four pictures displayed in Figure 6, that are obtained with the OQMC method using the same set of parameters, we can hardly detect any disparity in each pair of the corresponding sub-figures.

V. CONCLUSIONS AND DISCUSSION
In this work, we have investigated spin textures of a 2D triangular lattice in the presence of AF HE and chiral DM interactions with a quantum Monte Carlo method optimized recently. It is observed that within an external magnetic field applied normal to the 2D plane, well periodic and symmetric AF helical, skyrmionic and vortical lattices can be generated at temperatures higher than those predicted by previous authors; 19 the AF SL states prevail in a broad T − H region in the phase diagram; the sizes of skyrmions and vortices change discontinuously at a few critical points of the applied magnetic field; each of these AF periodical SLs and VLs can be decomposed into three FM sublattices; and the topological charge density of each FM sublattices of an AF SL also forms crystal which coincides with that FM SL almost exactly; the intensity of spin structure factor of each AF lattice varies in a very broad range in the 2D reciprocal space, where the bright points in the q x ∼ q y plane form symmetric patterns, but most points have very low intensity, which are connected continuously to form a blue background.
In fact, we have performed simulations with different D/|J | ratios, but only the results for D/|J | = 1 are presented in this manuscript for compactness. In general, the larger the ratio, the smaller periodicity of the AF SL; and our calculated T V L and T SL are usually several times larger than those estimated by the CMC method. 19 For instance, when D/J = -0. 5  To assess the results obtained with the OQMC mehtod, we have repeated the simulations by means of another quantum approach, namely the SCA method since the self-consistent algorithm is employed there. [21][22][23][24]26,27,29 The results obtained with J = -1, D = 1, and H = 7 are displayed in Figure S1 and S2 in comparison with those depicted in Figure 3 and 6, correspondingly.
The xy-projections of the AF SLs obtained at different temperatures using the two methods, when plotted together, overlap with each other almost exactly, as shown in Figure S1 In Figure S2 and Figure 6, the curves of the total topological charges versus changing temperature, the charge density distributions of the A, B and C FM sublattices of the AF SL at T = 0.05, that are calculated with the two methods, show no disparity.

Appendix B: Néel-Type AF SL Simulated with OQMC Method
By assigning J = -1, D = 1, and H = 7 which is perpendicular to the 2D antiferromagnet, OQMC simulations are started from a paramagnetic state, then carried out down to a very low temperature, T = 0.05. Now, D is assumed to be in-plane and normal to r ij , so the induced AF SL shown in Figure S3 and S4 are of Néel-type.
As shown in these two figures, this AF SL can also be decomposed to three FM SL, each of them contains 12 FM skyrmions in the 30 × 30 lattice, and the topological charge density of each FM sublattice forms a periodical crystal of the same pattern as the spin sublattice.
In comparison with the results shown in Figure 3 and 6, we can see that the two types of AF SLs obtained with the same set of parameters exhibit very similar characteristics. So the AF SLs of the Bloch and Néel types form a dual pair.