Vortex Structure in Magnetic Nanodots: Dipolar Interaction, Mobile Spin Model, Phase Transition and Melting
VVortex Structure in Magnetic Nanodots: Dipolar Interaction, Mobile Spin Model,Phase Transition and Melting
Aur´elien Bailly-Reyre and H. T. Diep LPNHE - Laboratoire de Physique Nucl´eaire et de Hautes ´Energies,Sorbonne Universit´e, Universit´e de Paris, CNRS/IN2P3,4 Place Jussieu, 75005 Paris, France ; [email protected] Laboratoire de Physique Th´eorique et Mod´elisation,CY Cergy Paris Universit´e, CNRS, UMR 8089, 2 Avenue Adolphe Chauvin,95302 Cergy-Pontoise, France; [email protected], corresponding author (Dated: February 9, 2021)We study in this article properties of a nanodot embedded in a support by Monte Carlo simulation.The nanodot is a piece of simple cubic lattice where each site is occupied by a mobile Heisenbergspin which can move from one lattice site to another under the effect of the temperature and itsinteraction with neighbors. We take into account a short-range exchange interaction between spinsand a long-range dipolar interaction. We show that the ground-state configuration is a vortex aroundthe dot central axis: the spins on the dot boundary lie in the xy plane but go out of plane witha net perpendicular magnetization at the dot center. Possible applications are discussed. Finite-temperature properties are studied. We show the characteristics of the surface melting and determinethe energy, the diffusion coefficient and the layer magnetizations as functions of temperature. PACS numbers: 75.10.-b ; 75.10.Hk ; 64.60.Cn :I. INTRODUCTION
In finite systems such as magnetic thin films and nan-odots, spin configurations are often non uniform near thesurface. The ground-state (GS) structure results fromthe competition between interactions in the system. Thecombination of the frustration [1] resulting from compet-ing interactions and the boundary effects in finite systemsgives rise to unexpected phenomena [2]. Among the com-peting forces, let us focus on the dipolar interaction whichfavors an in-plane non uniform spin configuration in flatand small samples, and the exchange interaction whichtends to align spins in parallel configuration. In addition,one can have the presence of a perpendicular anisotropyis known to arise with a large magnitude in ultrathinfilms [3, 4]. Note that, in thin films with Heisenberg andPotts models, the competition between the dipolar in-teraction, the exchange interaction and the perpendic-ular anisotropy causes a spin re-orientation transitionat a finite temperature [5, 6]. There has been a greatnumber of other works treating the dipolar interactionin the the presence of a perpendicular anisotropy in 2Dmonolayers and thin films. All of them found variousground states (GS) such as in-plane, out-of-plane andnon-uniform strip-domain configurations. Let us men-tion a few of them: in Ref. 7 an analytical calculationhas been performed at zero temperature ( T ) to find GSby varying the uniaxial surface anisotropy in a mono-layer and in thin films (i. e. infinite lateral dimension).In Refs. 8–10, Monte Carlo (MC) simulations have beencarried out for monolayers and thin films at finite T wherere-orientation phase transitions have been found. Ex-cept in Ref. 10 where all transitions are of second order,the two other works found first- and second-order transi- tions depending of the ratio of perpendicular anisotropyto dipolar strength. In Ref. 11, the authors used micro-magnetic simulations to calculate the T = 0 configura-tions for wires and disks. They did not consider finite- T behaviors. In Ref. [12], a spin-reorientation has beenexperimentally observed in Nd Fe B by using Lorentztransmission electron microscopy at variable tempera-tures and magnetic fields. It was shown that skyrmionsare created around the spin-reorientation temperature.The absence of works dealing with ultrafine dots at finite T using a mobile spin model has motivated the presentwork.In this paper, we focus on the case of a small magneticnanodot with mobile Heisenberg spins. Various GS con-figurations have been observed in such nanosystems [13],depending on the size of the sample, the ratio betweenthe exchange and dipolar interactions, and the type ofthe lattice. Systems in which the core vortex structureoccurs hold much promise from the commercial point ofview; the occurrence of this structure has already beendemonstrated experimentally [14–16] by different imag-ing techniques. A major advantage of core vortex struc-tures is the central region (core) of nonzero perpendicularmagnetization, the polarization of which is stable at roomtemperature as shown by Shinjo et al. [14]. Interestingly,core magnetization reversal [17, 18] can be realized in twoways, by applying a strong magnetic field perpendicular to the surface of the sample, or a short pulse of magneticfield parallel to it. This property of magnetic nanodotsopens the door to their application in magnetoresistiverandom access memory (MRAM).The current development of a technology that allowsto obtain samples with a very small dimension [19–21]has inspired us to investigate, with the use of MC simu-lations [22, 23], the behavior of the core vortex structure, a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b so interesting from the point of view of applications, un-der the impact of the temperature T using a mobile spinmodel. The mobility of the spins with increasing T givesthe opportunity to investigate the melting of these nan-odots.The GS structure found for a nanodot can be consid-ered as a single skyrmion [24]. There are several mech-anisms and interactions leading to the appearance ofskyrmions in various kinds of matter. The most popularone is certainly the Dzyaloshinskii-Moriya (DM) inter-action which was initially proposed to explain the weakferromagnetism observed in antiferromagnetic Mn com-pounds. The phenomenological Landau-Ginzburg modelintroduced by I. Dzyaloshinskii [25] was microscopicallyderived by T. Moriya [26]. The DM interaction has beenshown to generate skyrmions in thin films [27–29] and inmagneto-ferroelectric superlattices [30–32].In this paper, we study a magnetic nanodot embeddedin an non-deformable recipient, using the mobile Heisen-berg spins. The mobile Potts model with short-rangeinteraction has been used to study the phase transitionand the sublimation in a solid [33]. Here we extend thismodel to the case of a mobile Heisenberg model whichincludes a long-range dipolar interaction. The dot canbe heated to high temperatures to melt. Experimentally,one can imagine periodic arrays of such a dot embeddedon a crystal plane and transport properties of itinerantspins across such a plane can be studied in the presence ofdots. Itinerant spins are scattered by magnetic dots anddesired transport properties can be obtained by modify-ing dot arrangements.The purpose of this work is (i) to investigate the GSconfiguration in magnetic nanodots taking into accountthe short-range exchange interaction and the long-rangedipolar interaction, (ii) to study the nature of the or-dering and the phase transition at finite temperatures insuch nanodots, (iii) melting behavior. The methods weemploy in this paper are MC simulations with severaltechniques.The paper is organized as follows. Section II is devotedto the determination of the GS, while section III showsMC results of finite-temperature behaviors. Concludingremarks are given in section IV. II. GROUND STATEA. Mobile spin model
Let us consider a recipient of size N L = L × L × L z withthe cubic lattice where L is the x and y linear dimensions,and L z the size in the z direction. We consider a number N s of mobile Heisenberg spins which is less than the totalnumber of lattice sites N × N × N z . Therefore, eachlattice site can be occupied by a Heisenberg spin or canbe empty. The concentration of spins in the recipient istherefore c = N s /N L .We consider the following Hamiltonian containing the exchange interaction between nearest neighbors (NN)and the dipolar interaction between spins without cut-off: H = − NN (cid:88) i,j J ij S i · S j − D all (cid:88) i,j (cid:34) S i · (cid:126)r ij )( S j · (cid:126)r ij ) r ij − S i · S j r ij (cid:35) , (1)where J ij denotes the exchange integral between two NNspins i and j , D is the dipolar coupling parameter, S i ( | S i | = 1 for all i ) is the spin at the i -th site, and r ij ( r ij = | r ij | ) is the position vector connecting the spinsat the i -th and j -th sites. The dipolar energy is calcu-lated from the expression included in the Hamiltonian (1)without any numerical approximations; in particular wedo not introduce the cut-off radius, since this has beenshown [10, 13] to affect quantitatively the calculation re-sults in a sensitive manner. B. Method of ground-state determination
In the following, we take J ij = J ⊥ = 1 between NNin the z direction ( i.e. between NN in adjacent planes)and J ij = J // = 4 between in-plane NN. As will be seenbelow, another choice will not change our results but itchanges the range of values of D to have a vortex GS.We start from a random spin configuration of N s mo-bile spins. By random, we mean spins occupy random lat-tice positions and have random orientations. This statecorresponds to a high- T disordered phase. We slowlycool the system from that state using the Metropolis al-gorithm [22, 23]. The spins move from site to site, and atlow T , they condense into a film with a number of succes-sive layers fully occupied. This is due to the fact that thein-plane interaction is much larger than the perpendic-ular interaction. However, due to the long-range natureof dipolar interaction, the spin orientations are still ina weak disorder showing no symmetric configuration ex-pected from the symmetry of the Hamiltonian. To get ridof this and to come down to T = 0, we need to minimizelocally the energy of each spin as follows. We considerhere a spin localized on the lattice site at T = 0. To findthe ground state (GS) of the system we minimize the en-ergy of each spin, one after another. This can be numeri-cally achieved as the following. At each spin, we calculatethe local-field components acting on it from its NN usingthe above equations. Next we align the spin in its lo-cal field, i.e. taking S xi = H xi / (cid:113) ( H xi ) + ( H yi ) + ( H zi ) etc. The denominator is the modulus of the local field.In doing so, the spin modulus is normalized to be 1. Asseen from Eq. (1), the energy of the spin S i is minimum.We take another spin and repeat the same procedure un-til all spins are visited. This achieves one iteration. Wehave to do a sufficient number of iterations until the sys-tem energy converges. The spins at T = 0 form a dot ofsize L × L × L s which verifies c = N s /N L = L s /L z .An example of GS are displayed in Fig. 1 in a recipientof 15 × ×
12. The concentration used is c = 25%. Wesee that the dot size at T = 0 is 15 × ×
3. Let us givesome comments on Fig. 1. In Fig. 1a the dot is viewed in3D space. We see three compact layers. Figure 1b showsthe first layer structure, Fig. 1c shows the projectionon the xy plane where we can see that the projectionof the spins near the dot center are not exactly parallel.Figure 1d shows the projection on the xz plane, namelya side view of the dot. One sees that the spins at the dotcenter go out of the xy plane. This is easily understood inorder these spins reduce the spin orientation constraint.In a zero field in the z direction, these central spins canpoint in the + or - z direction. There is therefore a two-fold degeneracy. It may have application in magneticrecording memory using a small magnetic field along the ± z direction to control the magnetization direction.Note that there is a range of D where we observe suchvortex structure for a given ( J ⊥ , J // ). When we changethe recipient size and the concentration, this range of D changes. What we mean by this sentence is that when weincrease the system size for example, in the dipolar terma spin interacts with a larger number of spins becausethere is no cutoff in the dipolar sum of Hamiltonian (1).Therefore, to obtain a vortex which is favored by thedipolar interaction, we just need a smaller D to over-come the energy of the linear configuration favored by agiven set of ( J ⊥ , J // ). Another example of concentrationis shown in Fig. 2. Other cases displayed in Figs. 3-6below show that when we decrease the concentration c we have to increase D and vice-versa in order to havea vortex structure. It is not our purpose to preciselydetermine the range of D in which there is a vortex con-figuration. We give here only the physical mechanismwhich determines the planar vortex phase. III. MELTINGA. Slow heating
In this section, we slowly heat the dot in a recipient.We take the recipient height much larger than the thick-ness of the dot. The empty space in the z direction allowsthe dot to melt into a liquid at high temperature ( T ).Since the spins are mobile, a spin can move from onelattice site to a nearest site. MC simulation is used toupdate the spin position and the spin orientation. Thespin position is updated whenever there are vacant sitesnext to it, at the same time with the orientation update,using the Metropolis algorithm.We start the simulation using the GS state configu-ration and we heat the system to a temperature T . Asusual, we discard a large number of MC steps to equi-librate the system before averaging physical quantitiesover a large number of MC steps. Physical quantities FIG. 1. Ground state of c = 25% of mobile spins in a recipientof 15 × ×
12. The dot at T = 0 has the size 15 × ×
3. Wehave used J ⊥ = 1, J // = 4, D = 1: (a) 3D view, (b) First-layerconfiguration, (c) Projection on the xy plane, (d) Projectionon the xz plane. See text for comments. which are calculated include the energy E , the magneti-zation, the heat capacity and the magnetic susceptibilityas functions of T for different concentrations c , differentsizes of the recipient. We have also calculated the dif-fusion coefficient, which is the sum of the mean square FIG. 2. Ground state of c = 30% of mobile spins in a recipientof 15 × ×
12. The dot at T = 0 has the size 15 × × J ⊥ = 1, J // = 4, and with adipolar interaction D = 0 .
7: (a) 3D view, (b) Projection onthe xy plane, (c) Projection on the xz plane to see the spin z -components. The GS exhibits a vortex around the centerof the dot. The spins lie in the xy plane at the border exceptaround the vicinity of the center where they have non-zero z components. of the distance made by each spin at each T . We havealso computed the mean value of the number of nearestneighbors as a function of T . We recall some definitions:The total energy E = (cid:104)H(cid:105) (2)The Edwards-Anderson order parameter Q EA is Q EA = 1 N s ( t a − t ) (cid:88) i | t a (cid:88) t = t S i ( t ) | (3)where (cid:104)·(cid:105) indicates the thermal average with N s being thetotal number of spins. Note that the Q EA is calculatedby taking the time average of each spin before averagingover all spins of the system. This order parameter isvery useful in the case of disordered systems such as spinglasses [34] or doped compounds: it expresses the degreeof freezing of spins independent of whether the systemhas a long-range order or not [35].To study the change of behavior, we will show the re-sults of the following situations for comparison:i) all spins are supposed to be localized on their site(this case corresponds to the solid crystal)ii) spins on one, two, three, ... layers are mobile (thiscase allows us to study the partial melting process)iii) all spins are mobile. This case corresponds to thefull melting to the liquid phase at high enough T . B. Concentration effects
As said in the previous section, when we change theconcentration, we need another value of D to obtain avortex configuration. For 3 layers ( c = 25%), we pick D = 1 in the range of the vortex phase, while for 4 layers( c = 30%), we pick D = 0 . D for each concentration. Varying D slightly around thechosen values does not change qualitatively the results.Now if we work with one fixed variable, say c , and westudy the system behavior as a function of D , we find thevortex phase in a range of D , the collinear configurationsat small D , and no planar vortex at very high D . Thelinear phase and the non-planar vortex phase are notinteresting with respect to applications using the spinreversal by a small magnetic field.Let us show in Fig.3 the results of the energy as afunction of T for various dot thickness, namely variousconcentrations c with the choice of J // = 4 and J ⊥ = 1.All spins are mobile. FIG. 3. Energy per spin U vs T . Effect of the concentration.The lattice size of the system is 15 × ×
12 with J // = 4 and J ⊥ = 1. The blue curve corresponds to a spin concentrationequal to c = 50% with a magnetic dipole-dipole interactionequal to D = 0 .
3. The green curve corresponds to c = 30%and D = 0 .
7. The red curve corresponds to c = 25% and D = 1. It is interesting to note that the three energy curveshave the same transition temperature T c , a kind of fixedpoint independent of the concentration c . The Edwards-Anderson order parameter is shown in Fig. 4a confirmsthat the transitions of these concentrations occur at thesame temperature within statistical errors. We show inFig. 4b the diffusion coefficient. We see here that spinsstart to move after the transition at T c from the vortexconfiguration to the disordered phase.The fact that the meting temperature does not changewith the concentration up to 50% means that spins haveenough empty space to evaporate. We may imagine theextreme situation where c is close to 1: the first evapo-rated spins in the little empty space prevent the followingspins to evaporate due to the lack of empty sites. Thisincreases the melting temperature. FIG. 4. (a) The Edwards-Anderson order parameter Q EA vs T for several concentrations, (b) The diffusion coefficient C D vs T . The lattice size of the system is 15 × ×
12 with J // = 4 and J ⊥ = 1. The blue curve corresponds to c = 50%and D = 0 .
3. The green curve corresponds to c = 30% and D = 0 .
7. The red curve corresponds to c = 25% and D = 1. C. Comparison with a localized-spin system
We show here the case of a system of localized spinsto compare with the mobile spins shown above. Figure5 shows the energy for three thicknesses of the dot. Un-like the full melting case shown above, the transition oc-curs at a different temperature for a different thickness.This is well-known in magnetism of thin solid films: thetransition temperature increases with increasing thick-ness [35, 36]. We come back to this point at the end ofthis subsection.In addition to the Edwards-Anderson order parame-ter, we can define an order parameter in the case of nolong-range GS ordering: if the GS is well defined by anumerical method, then we can project the actual spinconfiguration of at a given T at a given time t on the GS.Needless to say, if the spin configuration is not stronglydeviated from the GS, the order parameter is close to 1.This is defines as: P ( T ) = 1 N s ( t a − t ) (cid:88) i | t a (cid:88) t = t S i ( T, t ) · S i ( T = 0) | (4)where S i ( T, t ) is the i -th spin at the time t , at tempera-ture T , and S i ( T = 0) is its state in the GS. The orderparameter P ( T ) is close to 1 at very low T where each FIG. 5. Energy vs T . Localized spin model. The latticesize of the system is 15 × ×
12 with J // = 4 and J ⊥ = 1.The blue curve corresponds to c = 50%, i.e. a 6-layer film,with a magnetic dipole-dipole interaction equal to D = 0 . c = 30% (4-layer film) and D = 0 .
7. The red curve corresponds to c = 25% (3-layer film)and D = 1. spin is only weakly deviated from its state in the GS. P ( T ) is zero when every spin strongly fluctuates in theparamagnetic state. The above definition of P ( T ) is sim-ilar to the Edward-Anderson order parameter used tomeasure the degree of freezing in spin glasses [34]: wefollow each spin with time evolving and take the spatialaverage at the end. However, the advantage of P ( T ) isthe fact that we can follow the GS configuration until itis broken.We show in Fig. 6 the Edwards-Anderson order pa-rameter Q EA and the projection order parameter P asfunctions of T . These confirm the difference of T c for dif-ferent film thicknesses. Note however that in magneticthin films with localized spins, T c increases with increas-ing thickness provided that all parameters other than thefilm thickness are the same [36]. In the present model,we cannot have this condition because we should choose D to produce the vortex structure while varying the filmthickness, explained previously in sections II B and III B.As a matter of fact, Figs. 5 and 6 show the data for threedifferent thicknesses with three different D , i. e. for threedifferent systems. The thickness is not the only variablehere. So the variation of T c is not due to it alone, thescaling in Ref. [36] does not apply. D. Effect on the number of layers that can melt
Let us show now the case where spins in one, two ormore layers are mobile. This allows us to follow the melt-ing progressively. Though artificial, this procedure cor-responds to a reality when the surface layer melts first,then the second layer, ... with increasing T [37].Figure 7 shows the energy and the diffusion coefficientversus T . The case of localized spins is also shown forcomparison. One sees that the transition temperaturedecreases as the number of mobile layers increases. Note FIG. 6. (a) Edwards-Anderson order parameter vs T , (b)Order parameter P . The model is the localized spin model.The lattice size of the system is 15 × ×
12 with J // = 4 and J ⊥ = 1. The blue curve corresponds to a 6-layer film, witha magnetic dipole-dipole interaction equal to D = 0 .
3. Thegreen curve corresponds to a 4-layer film and D = 0 .
7. Thered curve corresponds to a 3-layer film and D = 1. that the case of a system of completely localized spinshas a very high T c .Figure 8 shows the Edwards-Anderson order parame-ter and the order parameter P . Several remarks are inorder: (i) for a system of completely localized spins or asystem of completely mobile spins, there is only one tran-sition for each concentration, (ii) when a number of layersare mobile, the system has a partially ordered state: themobile layers become disordered at some T but the lo-calized layers are still ordered. This causes a step in theorder parameters observed in Fig. 8.Finally, we show in Fig. 9 the occupation rate R ofeach layer in the system in the case c = 30%. In the GS,and at low T ( <
1) only the first four layers are occupied( R R T , R of the first four layers diminish and R of the other layers increase. For very high T , the spinsoccupy all layers, and it is only after the melting that alllayers have the same occupation rate as expected in theliquid phase.To close this section, let us discuss about the values ofinteraction used above. The value J ⊥ = 1 was used as theenergy unit. The value J // = 4 was used to favor layeredordering. As long as J // > J ⊥ our results shown abovedo not change qualitatively, only the range of values of D giving rise to the vortex structure as well as the tran-sition temperature change. Finally, let us note that the FIG. 7. (a) Energy U vs T , (b) Diffusion coefficient C D vs T . The lattice size of the system is 15 × ×
12 with J // = 4, J ⊥ = 1 and D = 0 .
7. The concentration is fixed to c = 30%, i.e. results have been shown for the same lateral lattice sizebut changing this size will not change qualitatively theresults except the change of the range of values of D giv-ing rise to the vortex structure because of the long-rangenature of the dipolar interaction. IV. CONCLUSION
In this paper, we have studied a dot where the lat-tice sites are occupied by mobile Heisenberg spins. Thedot is embedded in a close recipient which allows to con-serve the number of spins. We have taken into accountthe in-plane and perpendicular exchange interactions andthe long-range dipolar interaction without cut-off. Theconfined geometry of the dot and the competition be-tween exchange interactions and the dipolar interactiongives rise to a ground state which is a vortex around thedot center with the spins at the center pointing out ofthe xy plane. Such a structure has a net perpendicularmagnetization with a two-fold degeneracy along the ± z axis. Using a small magnetic field one can pin the mag-netization in + or - z direction. If the dot is sandwichedbetween two ferromagnetic films, then one can obtain agiant magneto-resistance [38, 39] in a perpendicular spintransport. Let us denote the up spin of the ferromag- FIG. 8. System with partially mobile layers. (a) Edwards-Anderson order parameter vs T , (b) Projection of the orderparameter P vs T . The lattice size of the system is 15 × × J // = 4, J ⊥ = 1, and D = 0 .
7. The concentration is fixedto c = 30%, i.e. R per layer vs T .The lattice size ofthe system is 15 × ×
12 with J // = 4, J ⊥ = 1, and D = 0 . c = 30%, i.e. netic film by ↑ F , and the up dot magnetization by ↑ D .According to the giant magneto-resistance geometry, theconfiguration ↑ F | ↑ D | ↑ F will let the electron up spinsgo across the system (high current) while ↑ F | ↓ D | ↑ F will block the electron up spins (small current). Theswitch between the two states can be realized with anapplied magnetic field. Thanks to the smallness of thedot magnetization, one just needs a small magnetic fieldwhich does not heat the system with the magnetizationreversal. This is certainly an advantage over the use of aferromagnetic layer instead of the dot.Let us discuss another possible application of thepresent system in the domain of computer memory de-vices. The present system has a two-level structure whichis stable at finite temperatures below the melting point(cf. Fig. 4): center spins can be up or down piloted by anextremely small magnetic field to reverse just a few spinsas said earlier. This can serve as a two-bit unit. An ap-plication device one can imagine is a, array of dots withhorizontal lines are the bit lines and vertical lines are theword lines, similarly to what proposed in Ref. [40] usingarray of dots of skyrmions. The present system repre-sents a considerable advantage: small dot size (just a fewdozen of spins) and an extremely small magnetic field tooperate. It can therefore increase the stocking capacityin memory devices for example.We have studied the melting of such a dot with increas-ing temperature and found that, among other results,within the studied sizes the melting to the liquid phasetakes place at the same temperature regardless of the sys-tem size. This is not the case of the order-disorder phasetransition in a solid film [35]. Note that our techniqueusing the mobile spin model can be used to study the be-havior of liquid crystals. We have recently succeeded toobtain the nematic and smectic structures while coolinga liquid to low T using an appropriate choice of Hamil-tonian [41, 42].Finally, note that in this work, we have simulated thesystem at a fixed concentration using the canonical de-scription so that only one phase transition is observed.However, if we use the grand-canonical method, we be-lieve that we will observe the coexisting phases at a given T when varying the concentration, as what has beenfound in Ref. [33] using the mean-field approximation(cf. Fig. 6 of that reference). The implementation ofthe grand-canonical Monte Carlo method however is verycomplicated. This is left for a future study. [1] H. T. Diep (ed.), Frustrated Spin Systems , 3rd Edition,World Scientific, Singapore (July 2020).[2] H. T. Diep, V. Bocchetti, Danh-Tai Hoang andV. T. Ngo,
Theory and Simulation of Mag-netic Materials: Physics at Phase Frontiers ,http://arxiv.org/abs/1309.4754 (2013).[3] A. Zangwill,
Physics at Surfaces , Cambridge UniversityPress, London (1988).[4] J.A.C. Bland and B. Heinrich (editors),
Ultrathin Mag-netic Structures , vol. I and II, Springer-Verlag, Berlin(1994).[5] C. Santamaria and H. T. Diep, J. Mag. Mag. Mater. ,23 (2000).[6] D.-T. Hoang, M. Kasperski, H. Puszkarski and H. T.Diep, J. Phys.: Cond. Matter , 056006 (2013).[7] Y. Yafet and E. M. Gyorgy, Phys. Rev. B , 9145 (1988).[8] A. Moschel and K. D. Usadel, Phys. Rev. B , 16111(1995).[9] A. B. MacIsaac, K. De’Bell and J. P. Whitehead, Phys.Rev. Lett. , 616 (1998).[10] E. Y. VEdmedenko, H. P. Oepen, A. Ghazali, J.-C. S.L´evy and J. Kirschner, Phys. Rev. Lett. , 5884 (2000).[11] A. Maziewski, V. Zablotskii and M. Kisielewski, Phys.Rev. B , 134415 (2006).[12] Y. Xiao, F. J. Morvan, A. N. He, M. K. Wang,H. B. Luo, R. B. Jiao, W. X. Xia, G. P. Zhao,and J. P. Liu, Appl. Phys. Lett. , 132402 (2020);https://doi.org/10.1063/5.0022270.[13] C. S. Rocha, P. Z. Coura, S. A. Leonel, R. A. Dias andB. V. Costa, J. Appl. Phys. , 053903 (2010).[14] T. Shinjo, T. Okuno, R. Hassdorf, K. Shigeto and T.Ono, Science , 5481, 930-932 (2000).[15] J. Raabe, R. Pulwey, R. Sattler, T. Schweinb¨ock, J.Zweck and D. Weiss, J. Appl. Phys. , 4437 (2000).[16] A. Wachowiak, J. Wiebe, M. Bode, O. Pietzsch, M. Mor-genstern and R. Wiesendanger, Science , 18, 577-580(2002).[17] Q. F. Xiao, J. Rudge, E. Girgis, J. Kolthammer, and B.C. Choi, Y. K. Hong and G. W. Donohoe, J. Appl. Phys. , 103904 (2007).[18] N. Kikuchu, S. Okamoto, O. Kitakami, Y. Shimada, S.G. Kim, Y. Otani and F. Fukamichi, J. Appl. Phys. ,6548 (2001).[19] F. J. A. den Broeder, D. Kuiper, A. P. van de Mosselaerand W. Hoving, Phys. Rev. Lett. , 2769 (1988).[20] H. Kurt, M. Venkatesan and J. M. D. Coey, J. Appl.Phys. , 073916 (2010).[21] Y. Hodumi, J. Shi and Y. Nakamura, Appl. Phys. Lett. , 212506 (2007).[22] D. P. Landau and K. Binder, A Guide to Monte CarloSimulations in Statistical Physics , Cambridge UniversityPress, London (2009).[23] S. Brooks, A. Gelman, G. L. Jones and Xiao-Li Meng,
Handbook of Markov Chain Monte Carlo , CRC Press(2011). [24] T. H. R. Skyrme, Proc. Roy. Soc. A
127 (1961) ; Aunified field theory of mesons and baryons, Nucl. Phys. , 556 (1962).[25] I. E. Dzyaloshinskii, Thermodynamical Theory of ’Weak”Ferromagnetism in Antiferromagnetic Substances, Sov.Phys. JETP , 1259 (1957).[26] T. Moriya, Anisotropic superexchange interaction andweak ferromagnetism, Phys. Rev. , 91 (1960).[27] Sahbi El Hog, Aur´lien Bailly-Reyre and H. T.Diep, Stability and Phase Transition of SkyrmionCrystals Generated by Dzyaloshinskii- MoriyaInteraction, J. Mag. Mag. Mater. , 32-38(2018); https://doi.org/10.1016/j.jmmm.2017.10.031(2017); https://hal.archives-ouvertes.fr/hal-01474152;arXiv:1702.06841[28] H. T. Diep, S. El Hog and A. Bailly-Reyre, SkyrmionCrystals: Dynamics and Phase Transition, AIP Advances , 055707 (2018); https://doi.org/10.1063/1.5006269[29] H. T. Diep and H. Koibuchi, Frustrated Magnetic ThinFilms: Spin Waves and Skyrmions, chapter 10 in Frus-trated Spin Systems , 3rd Edition, Ed. H. T. Diep, WorldScientific, Singapore (July 2020).[30] I. F. Sharafullin, M. Kh. Kharrasov and H. T. Diep, Phys.Rev. B , 214420 (2019).[31] I. F. Sharafullin and H. T. Diep, Entropy , 862 (2020);doi:10.3390/e22080862.[32] I. F. Sharafullin and H. T. Diep, Symmetry , 26 (2020);doi:10.3390/sym12010026.[33] A. Bailly-Reyre, H. T. Diep, and M. Kaufman, Phys. RevE , 042160 (2015).[34] M. M´ezard, M. Parisi and M. Virasoro, Spin Glass The-ory and Beyond An Introduction to the Replica Methodand Its Applications , World Scientific (Singapore, 1986).[35] H. T. Diep,
Theory of Magnetism: Application to SurfacePhysics , World Scientific, New Jersey (2014).[36] X.T. Pham Phu, V. Thanh Ngo and H.T. Diep, Criti-cal behavior of magnetic thin films, Surface Science ,109-116 (2009), p. 114. doi:10.1016/j.susc.2008.10.037.[37] Virgile Bocchetti and H. T. Diep, Monte Carlo Simula-tion of Melting and Lattice Relaxation of the (111) Sur-face of Silver, Surface Science , 46 (2013).[38] M. N. Baibich, J. M. Broto, A. Fert, F. Nguyen Van Dau,F. Petroff, P. Etienne, G. Creuzet, A. Friederich, and J.Chazelas Phys. Rev. Lett. , 2472 (1988).[39] G. Binasch, P. Gr¨ u nberg, F. Saurenbach, and W. ZinnPhys. Rev. B , 4828(R) (1989).[40] Yadong Wang, Lei Wang, Jing Xia, Zhengxun Lai, GuoTian, Xichao Zhang, Zhipeng Hou,Xingsen Gao, WenboMi, Chun Feng, Min Zeng, Guofu Zhou, Guanghua Yu,Guangheng Wu,Yan Zhou, Wenhong Wang, Xi-xiangZhang and Junming Liu, Nature Communications ,3577 (2020).[41] A. Bailly-Reyre and H. T. Diep, Symmetry , 1574(2020); doi:10.3390/sym12091574.[42] V. Thanh Ngo, Phuong-Thuy Nguyen and H. T. Diep,Entropy22