Yu-Jiao Bo,Wen-Wen Li,Yu-Chen Guo and Ji-Chong Yang

Department of Physics,Liaoning Normal University,Dalian 116029,China

Abstract Compared with the traditional magnetic bubble,a skyrmion has a smaller size,and better stability and therefore is considered as a very promising candidate for future memory devices.When skyrmions are manipulated,erased and created,the density of skyrmions can be varied,however the relationship between the radii and the densities of skyrmions needs more exploration.In this paper,we study this problem both theoretically and by using the lattice simulation.The average radius of skyrmions as a function of material parameters,the strength of the external magnetic field and the density of skyrmions is obtained and verified.With this explicit function,the skyrmion radius can be easily predicted,which is helpful for the future study of skyrmion memory devices.

Keywords:radius of a skyrmion,shape of a skyrmion,lattice simulation

1.Introduction

Skyrmion is a topological soliton originally proposed to describe the baryons[1].In condensed matter,a particle-like object known as magnetic skyrmion was introduced theoretically in 1989[2].It was observed for the first time in 2D magnetic systems[3–6]involving Dzyaloshinskii-Moriya interactions(DMI)[7,8].Compared with the traditional magnetic bubble,the skyrmion is smaller,more stable and needs lower power to manipulate,therefore,it has been proposed that the skyrmion is a promising candidate for high density,high stability,high speed,high storage and low energy consumption memory devices[9–11].As a result,the magnetic skyrmions have drawn a lot of attention and been studied intensively recently[11–15].

A prerequisite for the use of skyrmions in devices is the knowledge of the relationship between the size of a skyrmion and parameters such as exchange strength,DMI strength and the strength of external magnetic field.Such a relationship can be investigated by solving the Euler–Lagrange equation of a skyrmion,for example numerically[16]or by using an ansatz[9],or by using the harmonic oscillation expansion[17],or by an asymptotic matching[18].It has been noticed that the radius of a skyrmion in the skyrmion phase is much smaller than that of an isolated skyrmion[17].Both the radii of an isolated skyrmion and the skyrmions in the skyrmion lattice were studied quantitatively in[19].In particular,numerical results were obtained for the equilibrium radii of skyrmion lattices.

However,as a potential candidate for storage,the skyrmion is meant to be manipulated,erased and created.In this case,the number of skyrmions can vary from only just one to filling the entire skyrmion lattice.The transformation of a skyrmion lattice to the saturated state is continuous,in this process,the skyrmion lattice gradually decomposes into isolated skyrmions in the saturated state[15].In this paper,we study the average radius of skyrmions with the density of skyrmions in the range between a single isolated skyrmion and the skyrmion lattice.While the results have been obtained for a single isolated skyrmion,and for skyrmion lattices,up to our knowledge,the radius of a skyrmion when the density of the skyrmions is between the skyrmion lattice and the single isolated skyrmion is poorly understood at a quantitative level.

The rest of the paper is organized as the following.The analytical and numerical results based on circular cell approximation are established in section 2.In section 3,we introduce the lattice simulation of Landau–Lifshitz-Gilbert(LLG)equation.We compare the theoretical results with the results of lattice simulation in section 4.A summary is made in section 5.

2.Circular cell approximation

The local magnetic moment of a skyrmion can be parameterized as

where r,φ,z are coordinates in a cylindrical coordinate,γ is the helicity angle,m=1 for a skyrmion and m=-1 for an anti-skyrmion,g=±1.The skyrmion number is Q=-mg.In the following,we only consider the skyrmion with Q=1(m=1,g=-1).

By using the circular cell approximation,the skyrmions are viewed as sitting in circular cells with radius R,which means the boundary condition θ(0)=π and θ(R)=0[19].In principle,θ(r)can be expanded using any Hilbert space.Since the wave-function of the ground state of the harmonic oscillator and the numerical solution of the Euler–Lagrange equation of a skyrmion are close in shape[17],we use the Hilbert space of harmonic oscillator to expand θ(r).We do not requireθ′(r)=0as in[17]becauseθ′(r)≠0is allowed by the Euler–Lagrange equation,therefore the eigen-functions of odd energy levels are also included.To impose the boundary conditions θ(0)=π and θ(R)=0,θ(r)to the next-to-next-to leading order can be written as

where φnare eigen-functions of harmonic oscillator,ω and c are parameters to be determined.By assuming the coefficients of the higher order terms are small,the power counting yields c～1/R2and 1/R≪1.

We concentrate on the case when the anisotropy is absent,the energy to be minimized iswith the energy density

where d≡D/J and b≡B/J,J is the strength of local ferromagnetic exchange,D is the strength of DMI,B is the strength of the external magnetic field which is assumed to be parallel to the z-axis.For simplicity,we consider dimensionless parameters,the matching is discussed in section 4.4.

Denoting s≡1/R,F can be expanded aswith

where γEis the Euler constant,Ci and Si are cosine and sine integral functions,andfm,n,nc,nsare constant numbers defined as

This seemingly lengthy expression ofnothing more than a polynomial ofand c.To achieve a higher precision,in principle,both the expansions of θ(r)and F can be worked out for higher orders.

For d=0.4,b=0.1,R=20,in the region that ω～0.15 and c～-0.01,we compare F within figure 1.can approximate F well in the region concerned.Especially,the positions where F andare minimized fit each other very well.To minimize F,we use variational method,so that there are two equations∂F/∂ω=0 and∂F/∂c=0,by which ω and c can be solved.

By setting a threshold h such that the sites with nz

Figure 1.Compare F withat d=0.4,b=0.1,R=20.The left panel is F andat ω=0.15 as functions of c,the right panel is F andat c=-0.012 as functions of ω.

It has been found that,by using the leading order ansatz θLO,to minimize F yields ω=w0b2/d2where w0=0.768548 is a constant[17].By using equation(6),the average radius rscan be written as

Figure 2.Compared with

For a skyrmion,nzvaries from-1 to 1 from the center to the edge,and the radius is determined by the number of sites with nz<1.But nzonly approaches 1,so we need to set a threshold slightly smaller than 1.If the approximate expression for θ(r)is sufficiently precise,the radius of the skyrmion is consistent as long as we choose a same h for theoretical predictions,lattice simulations and experimental measurements,and the verification of the theory is independent of the specific value of h taken.In this paper we choose 0.9 which is close to 1,other choices would lead to small differences,but would not change the conclusion.Then≈2.247 44,usingto approximate

In[19],the equilibrium R is numerically solved by minimize the energy density.Note that R in this case is independent of the density of the skyrmions.In our case,R is a quantity between the case of skyrmion lattice and the case of a single isolated skyrmion,and is determined by the density of the skyrmions,one can calculate rsafter R is given.

3.Lattice simulation

The lattice simulation is based on the LLG equation,denoting nras the local magnetic momentum at site r,the LLG can be written as[20–23]

where nris the local magnetic moment,α is the Gilbert damping constant and the effective magnetic field Beffis

with the discretized version of Hamiltonian defined as[24,25]

where δirefers to each neighbor.On a square lattice,one has δi=ei,therefore

The simulation was carried out on the GPU[20]which has a great advantage over CPUs because of the ability of parallel computing of the GPU.Equation(1)is numerically integrated by using the fourth-order Runge-Kutta method.

4.Numerical results

We run the simulation on a 512×512 square lattice.In the simulation,we use dimensionless homogeneous J,D and B.J=1 is used as the definition of the energy unit[25–27],the results are presented with d and b.In the previous works,the Gilbert constant was chose to be α=0.01 to 1[21,23,26–35].In this work,we use α=0.04 which is in the region of commonly used α.The time step is denoted asΔt.We use Δt=0.01 time unit,and the configurations typically become stable after about 106-107steps starting with a randomized initial state.The average radius of the skyrmions is measured as rs=As/NA where Asis the total area of the skyrmions which is determined by the number of sites in the isoheight nz=h with h=0.9,A=5122and N is the number of skyrmions.The standard errors of the radii of skyrmions are also measured.

To investigate the relationship between r,d and b,we simulate with d in the range of 0.2-0.6,and with growing b for each fixed d.We focus on those configurations that are in the skyrmion phase when stable.The phase diagram is shown in figure 3.

Figure 3.The phase diagram obtained by lattice simulation of LLG with randomized initial states.

4.1.The relationship between the average radius and densit y

In this subsection,we use the configurations at d=0.2,b=0.025,d=0.35,b=0.1,d=0.4,b=0.1,d=0.4,b=0.2,d=0.45,b=0.15 and d=0.6,b=0.15 to investigate the relationship between the average radius and density.Firstly we calculate the number of skyrmions in each configuration.By repeatedly and randomly erasing about 10%of the total number of skyrmions at a time and performing the simulation sequentially,we obtain the configurations at different N.Taking the case of d=0.4,b=0.1 as an example,the resulting configurations are shown in figure 4.

Because the size of the lattice is 512×512,The ratios of the average radii of skyrmions at different N to the radius of an isolated skyrmion are measured and denoted as.We compareandin figure 5.It can be seen that the theocratical resultsandcan approximately predictcorrectly,the deviations betweenandare generally within about 10%.Besides,for larger skyrmions,the theocratical results are generally better.Especially,for d=0.4,b=0.1,can fitvery well.For the cases where the sizes of skyrmions are relatively smaller,there are several possible reasons for the deviation.On one hand,when the skyrmions are smaller,they are not homogenously aligned,the distances between the skyrmions become larger,consequently the actual s is smaller for,so the points ofare biased towards a larger s.Secondly,the lattice simulation is more coarse for smaller skyrmions,which can also lead to differences with the theoretical results.Similarly,if the skyrmions are small and occupy only hundreds of sites,the θ(r)can no longer be treated as a continuous function.

Figure 4.The configurations corresponding to different N.

Figure 5.Comparewithand

Figure 6.Comparewithfor d=0.6,b=0.15.

Figure 7.Comparewithfor d=0.4,b=0.2.

Figure 8.rs at different d and b(marked as‘+’)and the fitted rs(d,b),i.e.Equation(10).

Figure 9.rs calculated with equation(9)(marked as‘+’)compared with equation(14)(the curved surface).

There are also cases that theare not fitted very well by the theocratical predictions.In the case of d=0.6,b=0.15,after erasing some skyrmions,the configuration began to enter the helical phase as shown in figure 6(a).In this case,when an isolated skyrmion is created,it is not stable and will grow into stripes.If we choose the average radius near the phase transition as a baseline,that is,we choose thenear the phase transition as 1,the results are shown in figure 6(b).It can be seen that theapproaches.

Another case is when d=0.4,b=0.2 as shown in figure 7.This is an example that when the skyrmions are rs=4.7003±0.0054,for N=103,rs=4.7005±0.0043,which are larger than the case of a signle isolated skyrmion rs=4.6865.One can see thathas a similar behavior.

4.2.A formula for the average radius in the skyrmion phase

Since we start from the randomized initial state,the obtained configurations are the skyrmion lattices.In this case,the relationship between rs,d and b are fitted by a rational function,the result is

small,they are no longer homogenously aligned.It can been found from figure 3 that the configuration is also near the phase transition from the skyrmion phase to the ferromagnetic phase.It is interesting that,the average radius is not always decreasing with s especially when s is small.For N=51,rsand fitted rs(d,b)(equation(10))are shown in figure 8.One can see that the rational function is consistent with the numerical results.

Figure 10.The shapes of isolated skyrmions in the saturated state.

We compare equation(9)with equation(14)in figure 9.Since the ρ for each configuration is different,the results of rs(b,d,ρ)are depicted as points in figure 9.Note that we remove the points correspond to the cases that the single isolated skyrmions are not stable such as d=0.6,b=0.15.It can be seen that rs(b,d,ρ)can match the results very well.rs(d,b,ρ)can be used to predict the average radius of the skyrmions when d,b and the density are given.

4.3.The shape of the skyrmion in the skyrmion phase

The function θ(r)is often used to describe the shape of a skyrmion[17].It has been assumed that the higher order corrections to θ(r)are small,which in fact requires that the shape of a skyrmion is not changed significantly in the skyrmion phase.Choosing the configurations at d=0.2,b=0.025,d=0.35,b=0.1,d=0.4,b=0.1 and d=0.45,b=0.15 as examples,we measure the average θ(r)with θ in the rangecos-1(0.9)≤θ≤cos-1(-0.9).The results are compared with θLO(r)with r rescaled according tofor example,in the case of d=0.4,b=0.1,one has s≈0.116 165 and θ(r,ρ)=θLO(r/0.783516).The results are shown in figure 10.

As shown in figure 10,the shapes of isolated skyrmions in a saturated state are similar to the shape of a single isolated skyrmion with r rescaled.This result indicates that our assumption is valid.Note that θ(r,ρ)is also the function θLO(r)with b rescaled asIt implies that the skyrmions can be seen as being experiencing an effective magnetic strengthwhen affected by other skyrmions.

4.4.Matching

In the calculations and lattice simulations,we use dimensionless parameters.The numerical results can be matched to the real material by using the rescaling introduced in[21,36].The rescaling factor is denoted as q andq=where λ is helical wavelength and a is the lattice spacing.The helical wavelengthes of real materials can be found in[12].For example,if we take λ≈60 nm,a=0.4 nm and D/J=0.4,then q≈6.75.Then r=6.2854 corresponds to r=6.2854×q×a≈16.97 nm.Meanwhile the time unit is rescaled aswhere J is the dimensionless exchange strength andJ′is the exchange strength of a real material.If we chooseJ′≈3 meV,the time unit ist′≈0.01 ns.The time step in the simulation isΔt=0.01t′≈1 ps.

5.Summary

One of the reasons that the skyrmion is proposed as a candidate for the future memory devices is because the size of a skyrmion is small.The radius of a skyrmion when the density of skyrmions is between the skyrmion lattice and the single isolated skyrmion is an important issue which is lack of exploration.In this paper,we study the average radius of skyrmions when the density is between a skyrmion lattice and a single isolated skyrmion.By using the harmonic oscillator expansion,the dependency of the average radius of skyrmions on the parameters of materials,strength of external magnetic field and the density of skyrmions is obtained theoretically.Then,a lattice simulation of LLG equation is performed to verify our results.

The theoretical result is presented in equation(9).Our result indicates that generally the average radius of skyrmions will decrease with the growth of density even when b and d are unchanged.The average radii at different d and b are measured by using lattice simulation.We confirm that our theoretical results can fit the simulated results well.With this relation,the skyrmion radius for different materials at different densities can be easily predicted.We also find that,the shapes of the skyrmions are insensitive to the density,which implies that the interactions between skyrmions can be seen as an effective magnetic strength.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China under Grant No.12 047 570 and the Natural Science Foundation of the Liaoning Scientific Committee Grant No.2019-BS-154.