Magnetic ordering and dynamics in monolayers and bilayers of chromium trihalides: atomistic simulations approach
S. Stagraczynski P. Balaz M. Jafari J. Barnas A. Dyrdal


Springs and a stopwatch: neural units with time-dependent multifunctionality
    Stephen Whitelam
    
============================================================================

empty



§ INTRODUCTION

Van der Waals magnetic (vdW) materials are currently of great interest due to expected applications in two-dimensional (2D) nanoelectronics and spintronics. These materials are built from weakly coupled 2D monolayers, and therefore they  can be relatively easily obtained in the multilayered form <cit.>, down to a single monolayer. The discovery of ferromagnetism in monolayers of CrI_3 <cit.> and   Cr_2Ge_2Te_6 <cit.> initiated an enormous interest in magnetic 2D vdW crystals <cit.>. There are several groups of well-known magnetic vdW materials, including chromium trihalides CrX_3 (X= I, Cl, Br, I) <cit.>, metal tribromides  MBr_3 (M=Mn, Cu, Fe, V), chromium-based ternary tellurides Cr_2XTe_6 (X = Ge, P) <cit.>, transition metal dichalcogenides MnX_2 and VX_2  (X = S, Se, Te) in two different phases denoted as 2H and 1T, <cit.>, and also others.  

The number of known 2D magnetic van der Waals materials is growing rapidly, and they have various structural, electronic, magnetic, and transport  properties. There is currently extensive search for new materials, with new properties, and new perspectives for applications, as these materials  seem to be ideal for building ultra-thin (of atomic scale) nanoelectronic devices (like spin valves or nonvolatile memory elements) for, e.g., quantum computing, spintronics, magnonics, opto-electronics, and others<cit.>. 
However, magnetic properties of most of the  known 2D van der Waals systems still require further investigations. The magnetic ground state of 2D vdW crystals depends, among others,  on their crystallographic phase, stacking geometry, substrate, and on possible twisting of adjacent monolayers. In many cases,  the magnetic ground state survives only at low temperatures, while at higher temperatures the systems become paramagnetic. Therefore, much efforts are focused on searching for materials with critical temperatures above the room temperature.    

Another group of 2D materials are layered molecular magnets. Indeed, a number of molecular magnets have been synthesized in the monolayer and bilayer forms <cit.>. Both 2D  vdW materials and 2D molecular magnets offer a unique opportunity to  test various theoretical models of magnetic interactions and magnetic physical phenomena in two dimensions. As for the interactions, this includes exchange symmetric coupling, Dzyaloshinskii-Moriya Interaction (DMI), magnetic anisotropy, higher-order interactions, like two-ion magnetic anisotropy (or Kitaeev term), and also others <cit.>. Many of these materials have the ground state spin orientation perpendicular to the monolayers or bilayers. In some cases, the magnetic anisotropy enforces in-plane spin orientation. Moreover, these materials also allow to observe various topological phases and phase transitions, like strictly 2D skyrmions or topological Berezinsky-Kosterlitz-Thouless phase transitions <cit.>.       

Among various groups of vdW materials, there are materials which are  promising for applications, like for instance transition metal dichalcogenides, especially Vanadium  dichalcogenides  VX_2 (X=S, Se, Te) <cit.>, and chromium  trihalides, e.g., chromium iodine and chromium chlorine, CrI_3 and CrCl_3
 <cit.>. In this paper we will focus on the latter group, specifically on CrI_3. Monolayers of this material are ferromagnetic, with the magnetic moments oriented perpendicularly to the layer plane <cit.>. This material is very interesting for  fundamental research <cit.>. However, the corresponding Curie temperature is around 40 K <cit.>, which  is too low for most of practical applications. In the bilayer form, the two monolayers are coupled antiferromagnetically. However  magnetic properties can be externally tuned by electric field <cit.> or strains <cit.>, and interlayer coupling can be even changed from antiferromagnetic to ferromagnetic one <cit.>. Moreover, magnetic properties of CrI_3 in a multilayer form depend on the stacking geometry <cit.>. This material seems to be ideal  for testing various theoretical models of spin configuration, including topological skyrmion textures <cit.>, as well as excited states (spin waves) which display certain topological features <cit.>. It was pointed out <cit.>, that DMI in CrI_3 is too small to generate skyrmion textures. However, it was also predicted that DMI in Janus monolayers of chromium trihalides Cr(I,X)_3 can be significantly larger than that in CrI_3, and can generate magnetic skyrmion textures <cit.>. Therefore, we assume a more general case, which admits larger values of DMI (e.g. in Janus structures).

In this paper we analyze magnetic properties of CrI_3, and to do this we employ simulation methods like atomistic spin dynamic (ASD) technique and Monte Carlo simulations <cit.>. From this we determine such magnetic properties like Curie temperature, hysteresis curves, (quasi)equilibrium spin configurations, or excited states (spin waves). We limit considerations to monolayers and bilayers of CrI_3. As some results on the monolayers are already available  in the literature, much less work has been done on the bilayer structures of  CrI_3. 
Atomic structure of the bilayers of CrI_3 is shown in Fig.<ref> for two different stacking of the monolayers. This figure shows top and side views. Atomic structure of a single  monolayer is effectively shown by the top part of Fig.<ref>(a).     



 


§ MODEL AND METHOD




 §.§ Structure and spin Hamiltonian


Magnetic properties of chromium trihalides remarkably depend  on the number of monolayers in the sample, and also on the geometry and offset of their stacking. This difference is especially significant for small numbers of layers in the system.  The most stable stacking geometries are the AA and Rhombohedral ones, see Fig.<ref>. The following analysis is limited to both monolayers and bilayers of CrI_3. Magnetic moments in CrI_3 are predominantly localized on the chromium atoms, while contributions from iodine atoms are rather small and may be omitted. However, the iodine atoms contribute to exchange interactions between the magnetic Cr atoms (e.g. via superexchange), and thus have a significant impact on magnetic properties of  CrI_3. The Cr atoms in an individual monolayer form a hexagonal (honeycomb) lattice, and positions of these atoms have been taken from the unit cell of bulk CrI_3.

Accordingly, magnetic properties of CrI_3 will be described by an effective spin Hamiltonian attributed to the corresponding hexagonal lattice of Cr atoms, which we split into two terms, 

    ℋ= ℋ_0 + ℋ_ DM,

where ℋ_0 includes the exchange interactions, magnetic anisotropies, and the Zeemann energy in an external magnetic field 𝐁.  The second term, ℋ_ DM, takes into account antisymmetric exchange, known also as Dzyaloshinskii-Moriya Interaction (DMI). 

Assuming isotropic  exchange interactions, we  write the first term in the form 

    ℋ_0 = -∑_ij J_ij𝐒_i ·𝐒_j - k∑_i (S^z_i)^2 - μ_S ∑_i𝐁·𝐒_i,

where the sum over i and j means the summation over lattice sites.  
Here, the positive  exchange parameters J_ij correspond to ferromagnetic (FM) coupling, whereas negative ones correspond to antiferromagnetic (AFM) coupling. In the following we take nonzero exchange parameters between nearest-neighbours, second-nearest and third-nearest neighbors. In a more general case the exchange parameters can be anisotropic. The second term describes the magnetic anisotropy. Here, positive  constant k corresponds to the  easy-axis anisotropy along the z-direction (normal to the layer), while negative k corresponds to the easy-(xy)plane anisotropy. Furthermore, μ_s is the atomic magnetic moment of Cr atoms. Here, and in the following, if not stated otherwise, the vector 𝐒 is a unit vector along the orientation of a local spin moment <cit.>. Accordingly, one needs to adjust  the exchange parameters, which usually are related to spin Hamiltonian, where 𝐒 is the spin operator.     

The second term in Eq.(<ref>) includes the  Dzyaloshinskii-Moriya interactions (or alternatively antisymmetric exchange). Existence of nonzero components of this interaction depends on symmetry of the system. In van der Waals materials one can easily tune externally the symmetry and thus also the components of DMI, e.g., by an external mechanical strain or by electric field normal to the layers. In the case under consideration, the DMI  Hamiltonian can be written in the form <cit.>  

    ℋ_ DM =    -∑_⟨ i, j ⟩𝐝^NN_ij·( 𝐒_i ×𝐒_j )          -∑_⟨⟨ i, j ⟩⟩𝐝^NNN_ij·( 𝐒_i ×𝐒_j ),

where the first term includes interaction between intralayer nearest-neighbours (NN), while the second term between next nearest-neighbours (NNN). Here, 𝐝^NN_ij and 𝐝^NNN_ij are the corresponding Dzyaloshinskii-Moriya vectors <cit.>.      
This vector for NN has the following form:  𝐝^NN_ij = d^NN_xyẑ×l̂_ij + τ_ij^NN d^NN_z ẑ, where ẑ is a unit vector along the axis z (normal to the layer), l̂_ij is the unit vector oriented from the site i to site j, while d^NN_xy and d^NN_z are parameters  and τ_ij^NN =1 (-1)  for the AB (BA) link. 
In turn, in the second term 
𝐝^NNN_ij = τ^NNN_ij( d_xy^NNNl̂_ij + ν_ij d_z^NNNẑ), with τ^NNN_ij = +1 for AA-link  and -1 for BB-link, while   ν^NNN_ij = ± 1 and alternates its sign on the consecutive NNN links of a chosen site <cit.>. 

The DMI between nearest-neighbors in pristine CrI_3 vanishes, while the DMI  between the next-nearest-neighbors is  rather small, d_z^NNN=-8.8μeV and d_xy^NNN=73.0μeV <cit.>.
However, DMI can be tuned by an external electric field normal to the layer, so the parameters d_xy and d_z can be nonzero and relatively large <cit.>. It can be also tuned by external strain. Moreover, the DMI can be  significantly enhanced in the corresponding  Janus monolayers of chromium trihalides Cr(I,X)_3 <cit.>.  All this makes the DMI parameters between NNs and NNNs tunable in a quite considerable range. Finally,we 
note, that in the Hamiltonian assumed here we neglect such  terms, like anisotropic exchange and higher-order anisotropy terms (e.g., two-ion anisotropy or Kitaev term).
     
    


 §.§ Numerical procedures




Our numerical calculations are based on the  atomistic spin dynamics (ASD) method implemented into the Vampire code <cit.>. This method allows for taking into account all details of the system geometry, and also for keeping insight into behaviour of the spin moments during the simulation process. The ASD simulations are carried out by solving the coupled system of atomistic stochastic (at finite temperature) Landau–Lifshitz-Gilbert (LLG) equations, 

    ∂𝐒_i/∂ t = -γ/1 + α^2[𝐒_i ×𝐁^eff_i + α𝐒_i× (𝐒_i ×𝐁^eff_i )],

where γ is the gyromagnetic ratio, α stands for the Gilbert damping parameter, and 𝐁^eff_i is the effective magnetic field acting on the spin at site i.
This effective field 𝐁^eff_i consists of two terms: one is deterministic and is related to the spin Hamiltonian, while the second one is a stochastic term, that is determined by thermal fluctuations, 

    𝐁^eff_i = -1/M_s∂ℋ/∂𝐒_i + ξ^ (th)_i.

In the low energy regime, the stochastic term, ξ^ (th)_i is usually modelled by the Gaussian thermal noise. 


To verify the results on spin dynamics obtained by the atomistic LLG equations, we also  used other methods,  like the method based on the spin-spin correlation function and  the  Monte Carlo simulation technique. We find a very good agreement between the results obtained by the atomistic LLG equations and those obtained by the Monte Carlo simulations. As the former technique is more appropriate for dynamical properties, the latter one is more convenient for static properties. 
However, the static properties are also accessible from spin dynamics, but require much higher computational efforts.  In addition, to study dynamical properties we also use the atomistic approach implemented within the Uppsala ASD code <cit.>.
 
 An important step in the simulations process, especially when simulating    
 hysteresis loops, critical temperatures, or magnon dispersion relations, is the preparation of  the relevant initial spin configuration. To do this we use the field-cooling procedure. 
 We  start this procedure always with the random spin state at a high temperature (above the Curie temperature). After performing one million of Monte Carlo equilibration steps (or time steps for LLG equation), the temperature becomes decreased nearly adiabatically towards the minimum temperature (i.e. temperature, at which the required simulations are to be performed, e.g. 15 K in most of our simulations of the hysteresis loops). 

 In the following we will apply the above mentioned methods to determine the Curie (Néel) temperature, as well as hysteresis loops at a finite (nonzero) temperature. The simulations will be performed for monolayers, as well as for bilayers in two different stacking geometries, i.e., in the AA and rhmbohedral configurations. In addition we will also consider dynamical properties of the systems under consideration, especially spin wave excitations.    
 


§ RESULTS
 



 §.§ Critical temperatures

The reported experimental data indicate on a ferromagnetic  ground state in CrI_3 monolayer <cit.>. Moreover, FM ordering  is also expected in multilayers of chromium trihalides, (i.e. in corresponding bulk material). However, the bilayer of CrI_3  was reported to display antiferromagnetic interlayer exchange coupling <cit.>, so such bilayers should also reveal not only features typical of ferromagnets, but also some properties typical of antiferromagnets.
The interlayer coupling can be tuned (including also sign change), e.g., by an external gate voltage <cit.> or strain<cit.>.  Therefore, one of our main objectives in the case of bilayers is to  investigate the impact of the layer stacking on the critical temperature. The analyzed geometries (AA and rhombohedral) are visualized in Fig. <ref>. 
Second, in antiferromagnetically coupled bilayers one may expect certain features related to the breakdown of ferromagnetic order within the monolayers and of antiferromagnetic order related to the interlayer coupling (Curie and Néel temperatures.
    

    
The simulated data (points) for both monolayer and ferromagnetic bilayer (ferromagnetic interlayer coupling) in the AA and Rh stacking geometries are shown in Fig. <ref>(a). We included there  the exchange parameters between nearest neighbours, J_1, second-nearest neighbours, J_2, third-nearest neighbours, J_3, and interlayer exchange coupling between nearest neighbours, J_z, as described in the figure caption. The other parameters are also described there.       
To determine the  Curie temperatures, the corresponding data have been fitted to the formula (1-T/T_C)^β, and from this fitting we found T_C ≈ 33.93 K for the monolayer, T_C ≈ 37.64 K for the Rh-stacked bilayer, and T_C ≈ 39.64 K for the AA-stacked bilayer. The Curie temperature for the monolayer is lower than that for the bilayers, and the bilayer in the AA-stacking has Curie temperature higher than that for the Rh-stacking. For the assumed parameters, the DMI has no significant impact on the magnetization vs. temperature profile and on Curie temperature. 

The critical temperatures (both Curie and Néel) also can be determined from magnetic susceptibility or specific heat. In the case of ferromagnetic monolayer and ferromagnetically coupled bilayers, both these methods give results which are consistent with those obtained from the above magnetization analysis. Therefore, to determine the critical temperature in bilayers with antiferromagnetic interlayer coupling, we  used the  Vampire package <cit.> to simulate the specific heat as a function of temperature, C_v = [⟨ U^2 ⟩ - ⟨ U ⟩^2]/k_B T^2, where U is the internal energy. The specific heat has a cusp at the critical temperature, and this property is used to determine the Curie temperature, T_c, in ferromagnetic systems, and the Néel  temperature in antiferromagnets. 

The numerical results on the specific heat for  bilayers with antiferromagnetic coupling in the AA and Rh stackings and also for the ferromagnetic monolayer 
are shown in Fig. <ref>(b). The critical temperatures  for the bilayers with antiferromagnetic interlayer coupling are clearly seen and are roughly the same as Curie temperatures for the corresponding bilayers with ferromagnetic interlayer coupling. This is because absolute value of interlayer coupling is remarkably smaller than the intralayer exchange coupling constant J_1, and absolute magnitude of the interlayer coupling is assumed the same in bilayers with ferro- and antiferromagnetic interlayer coupling. The curve for ferromagnetic monolayer is presented in Fig. <ref>(b) just to show that the methods used in simulations of the results in Fig. <ref>(a) and Fig. <ref>(b)  give the same critical temperatures. 

When the interlayer exchange coupling is antiferromagnetic, the total magnetic moment of the bilayer vanishes in the ground state. As already mentioned above, the absolute magnitude of the interlayer coupling is remarkably  smaller than the intralayer exchange coupling J_1. Thus, with increasing temperature, the antiferromagnetic moment (Néel vector) vanishes first and we have a system of two ferromagnetic layers (with average moments either aligned or antialigned). When the temperature grows further the magnetic moment in each monolayer (and thus in the whole system) drops to zero at the Curie temperature as determined above. Thus, on the simulated magnetization vs temperature curves one could expect two cusps, one at a lower temperature that corresponds to the breakdown of the antiferromagnetic alignment of the two monolayers, and another one at higher Curie temperature. However, in the simulation data on the magnetization vs temperature, there is no signature of the breakdown of antiferromagnetic configuration. Apparently, it would be resolved in a multilayer system <cit.>. 



In order to get some information on the temperature at which the antiferromagnetic alignment of the two ferromagnetic monolayers is destabilized by thermal fluctuations, we analyzed average values of the total spin z-component of each monolayer separately as a function of temperature.

Results of the simulations are shown in Fig. <ref>, where the simulated normalised values of S^z for the two monolayers are presented  with blue and yellow dots, while the average values are shown by the solid blue and red lines.  With increasing temperature, fluctuations around the average values become larger, and at a certain temperature (recognized as the relevant Néel temperature) the average spin moments of the two monolayers start to change sign. The Néel temperature was identified as a temperature at which spin moments of the monolayers begin to fluctuate between the up and down orientations, see Fig. <ref>. The corresponding  Néel temperature is lower than the Curie temperature by about 3.5 K, as one might expect.    



 §.§ Hysteresis loops and spin textures


    
Hysteresis curves are important characteristics as they include information on magnetization variation with increasing external magnetic field and thus also on possible field-induced phase transitions. 
We have determined the hysteresis loops from the Monte-Carlo (MC) simulations as well as from the atomistic Landau-Lifshitz-Gilbert (LLG) equation. The starting point for each simulation was a random state at temperature T=80K and cooling field H_fc=1T applied perpendicularly to the sample plane. The field-cooling procedure was used to cool the system's temperature down to T=15K, and then the simulations of a hysteresis loop was initiated.
 


  §.§.§ Monolayers

    

 
Figure <ref> shows the hysteresis curves for fixed values of DMI between NNNs (as described above in the text), and different values of the  components d_z and d_xy of the DM vector between NNs. When d_xy=0 , then increasing the component d_z leads mainly to reduced  coercivity of the hysteresis curves, while their shape remains similar to those for ferromagnetic layers, see Figure <ref>(a). 
A more interesting impact of DMI on the hysteresis curves appears for the d_xy component of the DMI vector. The corresponding hysteresis curves for monolayers are shown in  Figure <ref>(b)  for different values of  d_xy and d_z=0. For small values of d_xy/J_1, the hysteresis curves are qualitatively similar to those for zero d_xy, but the corresponding coercivity becomes reduced with increasing d_xy. When d_xy grows further, some new features appear in the curves, see e.g. the curves for d_xy/J_1=0.3 and d_xy/J_1=0.35, which clearly distinguish these loops from those typical for ferromagnetic systems.  
The interplay of both components of NN DMI vector is presented in Figure <ref>(c), 
where both components are nonzero and comparable. 
For the assumed parameters the  corresponding hysteresis loops  contain clearly  features of both  d_z  and d_xy components of the DMI vector.    
From the hysteresis curves shown above one may conclude that DMI modifies magnetic texture of the system and this modification leads to the observed features in the hysteresis loops. In the curves shown here, the most pronounced features are for d_xy =0.35J_1. Therefore we analyse now in more details this specific situation (d_z=0 and d_xy =0.35J_1), and for clarity, we also omit here the NNN DMI.  
   
    
    
        

  
    In Figure <ref>(a) we show the average hysteresis curve  (solid line) and  the standard deviation from the average (shaded areas)  
    for d_xy =0.35J_1, d_z =0 and 15 K.  
    The corresponding spin patterns for the points A, B, C, D, E and F, indicated on the hysteresis curve are shown in Figure <ref>(b). 
    

    Here, the results of MC simulations are shown in the top row. Then, the final set of MC simulations was used as initial step for simulations based on the atomistic LLG equation. The middle row shows the magnetic textures after 1ns of spin-dynamics using the LLG Heun method with damping factor λ=0.1. We have changed the numerical integrator from Monte Carlo to LLG in order to check if the spin texture was stable under both integrators. 
    
    The results from the middle row were subsequently used in the next simulation for the cooling temperature from 15K down to  0K.
    As the patterns for T=15K  have spin textures irregular at the boundaries due to thermal fluctuations, the edges become regular when cooling the system temperature down to 0 K, see  bottom row in Figure <ref>(b). As one can see, upon cooling the edges of various objects become smooth and one can clearly distinguish there regular circular objects and stripe structures. The circular objects are Néel-type  skyrmions. To see this in more details, in  Figure <ref>(c) we zoomed out an individual circular object, so now the spin structure associated with  the skyrmion is  clearly visible. To show further details, on the right side of this figure we present a crossection through the center of the  skyrmion . 




  §.§.§ Bilayers


Hysteresis loops in bilayers  depend  of the nature of interlayer coupling, and for ferromagnetic and antiferomagnetic couplings are usually remarkably different. In addition, they also depend on the stacking geometry. The hysteresis curves for the AA and Rhombohedral stackings in the absence of DM interactions and for isotropic exchange interactions are shown in Fig.<ref>(a) for  systems with ferromagnetic interlayer coupling and in Fig.<ref>(b).
for systems with antiferromagnetic interlayer coupling. 

    
    

    From this figure follows that the hysteresis curves for ferromagnetically coupled bilayers are similar to  typical hysteresis loops  of ferromagnetic layers,  and are roughly independent of the stacking geometry. On the other hand, the hysteresis curves for bilayers with  AFM interlayer coupling are significantly different from those for ferromagnetically coupled bilayers, and also different for different stacking geometries.   
    An important feature of the hysteresis loops for antiferromagnetic interlayer coupling is a wide plateau  at zero magnetization, which corresponds to antialigned magnetizations of the two monolayers. At a certain value of magnetic field there is a jump in magnetization associated with transition to the spin-flop phase (note, the magnetic field is normal to the layers). Then the system goes smoothly to the saturated state. The magnetic field, at which the transition to the spin-flop phase occurs with increasing magnetic field is different from the field at which the system goes back to antialigned configuration when the magnetic field  decreases. This leads to minor loops clearly visible in Fig.<ref>(b) on both sides (positive and negative magnetic fields). 
   These loops follow from the magnetic anisotropy in the system, and also depend on the stacking geometry. No such features appear in ferromagnetically coupled bilayers. The saturation field also depends on the stacking configuration. This is because the saturation field depends on the effective coupling between the two monolayers, and this coupling depends, in turn, on the  strengths of interlayer exchange parameter, and on the number of interlayer nearest neighbours. Note, the hysteresis curves for the bilayers in the AA and Rh stacking have been calculated for the same interlayer exchange parameters. However, the numbers of nearest neighbors in the Rh stacking is larger (but distance between them is longer then in the AA stacking). Thus, for the assumed parameters the  saturation fields are different for  both stacking geometries and is larger for the Rh stacking.             

   
 The impact of both NN and NNN terms of DMI  on the hysteresis curves in the Rh stacking geometry and for ferromagnetic and antiferromagnetic interlayer coupling is shown in Fig.<ref>(a) and Fig.<ref>(a), respectively, while the corresponding spin textures are shown in Fig.<ref>(b) and Fig.<ref>(b). As already discussed above, the DMI between NNNs is rather small in pristine CrI_3. In turn, the DMI between NN's can be easily tuned by external gate and therefore can be dominant (even if it vanishes in pristine CrI_3), so in Fig.<ref>(a) and Fig.<ref>(a) we assumed intrinsic values of  NNN DMI (as given in the text) and quite significant value of DMI between NNs, d_xy = 0.35 J_1 and d_z=0. We focus on the impact of the component d_xy of DMI as this component is crucial in the skyrmion texture formation.           
As in monolayers, the DMI modifies the hysteresis loops and also the corresponding spin patterns. These patterns include areas of various shapes, with irregular edges due to thermal fluctuations at finite temperatures.   However, when the temperature is cooled down to T=0K, these edges become regular and reveal among others stripe structures and well defined skyrmion textures. 
 



    
  
Consider first the bilayers with ferromagnetic interlayer coupling. The corresponding spin patterns are shown in Fig.<ref>(b)  for  a few points indicated on the corresponding hysteresis curve, see Fig.<ref>(a). In general, we find there stripes, skyrmions, and also other structures. Only spin patterns of the top monolayer are shown there as the spin patterns of the bottom monolayer are practically the same due to ferromagnetic interlayer coupling. Certain, almost not resolved difference may appear due to a shift of the two monolayers in the Rh stacking. Accordingly, the  skyrmions in both layers are coupled ferromagnetically and have the same topological charge. 


In bilayers with antiferromagnetic coupling, the corresponding spin patterns are qualitatively similar to those in FM-coupled bilayers, and we find the spin stripes, wavy stripes, and skyrmions, as shown in  Fig.<ref>(b). Though they look qualitatively similar to those in FM-coupled bilayers, they are different and this difference originates from different magnetic fields for which the spin patterns have been simulated and also from opposite signs of the interlayer coupling.  
IN the simulation results shown in Fig.<ref>(b), we show the spin patterns of both top and bottom layers. Note, for small magnetic fields (H=0 and H=0.5 T), the spin patterns of  both layers are practically opposite due to antiferromagnetic interlayer coupling. The situation is different for stronger magnetic field, H=2.28 T, which is strong enough to overcome te interlayer exchange coupling. In that case  the spin patterns of both layer reveal individual skyrmion in the spin pattern of the top layer and also individual skyrmion in the bottom layer. These two skyrmions, however, are not coupled and appear in different positions.  

Finally, we
also note  here, that the hysteresis loops and the corresponding spin patterns of the monolayers and bilayers, as shown above,  are for isotropic exchange interactions. These loops and spin patterns may become remarkably modified when the exchange interaction has anisotropic contributions, or more complex anisotropy terms (like two-ion anisotrop) occur in the system. A typical example  of such interactions is the so-called Kitaev term.  



 §.§ Artificially created Skyrmion lattices
 

In the spin patterns shown above for monolayers, Figure <ref>, and for bilayers,   Fig.<ref>(b) and Fig.<ref>(b), only single skyrmion states have been found in the corresponding simulation processes for the assumed parameters. Of course, for other parameters, especially for stronger DMI, one can achieve stable multi-skyrmion textures or skyrmion lattices <cit.>. However, we show now that one can also create arrays of (quasi)stable skyrmions by appropriate artificial arrangement of initial spin states, by reversing manually the local magnetization at the nodes of a certain lattice, see the top panel in Fig.<ref>.  The parts (a) and (e) of this panel present random initial states, while the others are arranged manually in the magnetic field 0.3T (b-d) and  0.748T (f-h). The evolution of the initial states is analyzed by Monte Carlo simulations, and the second and third rows of these columns show the spin structure after  300k and 1M of Monte Carlo steps. 

    

In case (a) we get at the and (after 1M steps) individual as well as coupled skyrmions. In the case (b) the initially created square lattice evolves into a square lattice of skyrmions with a size smaller than the size of the initially created areas of reversed magnetization. In turn, in case (c) and (d) the system after 1M of simulation steps evolves to smaller skyrmions arranged in a defected hexagonal lattice.
The origin of the defects might be twofold. First, they appear naturaly, because of the square shape of the simulation cell, which is incompatible with a hexagonal lattice. Second, defects might be excited at increasing temperature.

The columns (e) to (h) in Fig.<ref> correspond to the same initial states as in columns (a) to (d), but simulations are in the magnetic field of 0.748T. In all the four cases (e) to (h), after 1M of Monte Carlo steps we get the arrays of skyrmions (of approximately equal size). Interestingly, the square lattice in column (f) survives after 1M of simulation steps. 
From the above results one may conclude that the simulated skyrmion lattices may be stable or quasi-stable in time, or can evolve into a disordered hexagonal lattice. This also shows that one can intentionally create skyrmion lattices, which may be of some interest in skyrmionics. However, their stabilization requires further investigations <cit.>  



 §.§ Dynamical properties
 

Now, we consider briefly the spin excitations (spin waves) in both monolayers and bilayers of CrI_3. Therefore, we assume zero magnetic filed to exclude skyrmin texture formation when DMI is relatively large. First, we consider spin waves in bilayers assuming intrinsic NNN DMI parameters as given in the main text. As for the NN DMI we assume rather small its component d_xy, i.e., d_xy = 0.05J_1, which may be induced externally, e.g, by a gate voltage or strain due to a substrate, while for the second component we assume d_z=0. We focus on spin waves in bilayers in the AA and Rh stackings and with FM and AFM interlayer exchange coupling. 

The corresponding dispersion curves for the  bilayers with FM and AFM interlayer exchange coupling are shown in   Fig.<ref>.  The dispersion curves have been calculated using ASD Vampire code <cit.>, and  are shown there  along the main paths in the Brillouin zone. 
Within this approach, dynamics of individual spins is described in terms of Landau-Lifshitz-Gilber equation. To determine the spin wave mode, the system is externally excited and then Fourier analysis of the time evolution is used to determine the spectrum. 

As in a monolayer there are two modes <cit.>, one can expect four modes in bilayers due to splitting induced by the exchange coupling between the two individual monolayers. However, the interlayer coupling is rather small, so the splitting is not well resolved. 
Indeed, in the case of FM interlayer coupling, the splitting of the modes are not resolved in  Figs<ref>(a,b) for the AA and Rh stackings. The splitting is clearly resolved in the case of AFM interlayer coupling in the AA stacking, see Fig.<ref>(c), while it is not well resolved in the Rh stacking, see Fig.<ref>(d) However, in the latter case there is a small gap in the spectrum, induced by the NNN DMI. However, this gap is very small  due to small values of NNN DMI parameters. To check variation of this gap with NNN DMI  we consider the spin wave spectrum in the monolayer case, where the NNN DMI  can be significantly enhanced, e.g., in the corresponding chromium trihalide Janus structures <cit.>.        





Now we analyse spin waves in  monolayers <cit.>, assuming  relatively large values of the NNN DMI parameters, and focus on the topological features in the spectrum  associated  with NNN DMI. Though this term in DMI is small in pristine CrI_3, it can be remarkably enhanced, e.g., in the chromium trihalide Janus structures <cit.>, to produce remarkable effects in the spin wave spectrum by creating a gap in the Dirac point K of the Brillouin zone. 
The  spin wave modes for monolayers are simulated within the ASD Uppsala code <cit.> . As for exchange parameters, we note, that these parameters are not well determined yet and there is some scatter in the available data <cit.>. We assume here the parameters consistent with those provided in Ref. <cit.>. 
The results of simulations are shown in Fig.<ref> for four different values of the NNN DMI parameters.
In the absence of NNN DMI, there is no gap in the spectrum. However, with increasing DMI parameters the gap opens in the spectrum at the Dirac point K of the Brillouin zone, and increases with increasing parameters of NNN DMI.  



§ SUMMARY AND CONCLUSIONS


In this paper we have studied magnetic properties of chromium trihalides using numerical simulations within the atomistic spin dynamics and Monte Carlo  methods. These simulation techniques have been used to determine  critical  temperatures, hysteresis curves, and spin wave frequencies. These characteristics have been analysed for monlayers as well as bilayers in two (AA and rhombohedral) stackings and for ferromagnetic and antiferromagnetic interlayer exchange coupling. The main focus was on the influence of Dzyaloshinskii-Moriya interaction on the basic magnetic properties, especially on the hysteresis curves, spin textures (stripe domains, skyrmion textures, and others), and on the spin wave excitations, where Dzyaloshinskii-Moriya coupling between next-nearest-neighbours leads to opening a gap in the spectrum at the Dirac K point. Both, skyrmion formation and gap opening require quite a relatively strong DMI between nearest-neighbours and next-nearest-neighbours. Accordingly to observe these features one needs enhance the DMI either by external electric field (or strain), or use chromium trihalide Janus structures <cit.>. 










§ ACKNOWLEDGEMENTS (NOT COMPULSORY)


This work has been supported by the Norwegian Fi- nancial Mechanism 2014–2021 under the Polish Nor- wegian Research Project NCN GRIEG 2Dtronics no. 2019/34/H/ST3/00515.



§ AUTHOR CONTRIBUTIONS STATEMENT


S.S. performed numerical simulations using Vampire code. P.B. performed numerical simulations using ASD Uppsala code. M.J. supported numerical calculations of S.S.. A.D. conceived of the idea. A.D and J.B supervised the work and wrote the manuscript with
input from all authors. 








