Skyrmion glass in a 2D Heisenberg ferromagnet with quenched disorder

A 2D Heisenberg ferromagnet with exchange J and random magnetic anisotropy of strength D R ≪ J has been studied. Analytical theory for the dependence of the average size of a pinned skyrmion on the magnetic field H, and for stability of such skyrmions on a lattice, has been developed. It has been complemented by numerical studies of 2D lattices containing up to 40 million spins. At low fields the average size of the skyrmion, λ, is determined by the average size of Imry–Ma domains. On increasing the field the skyrmions first shrink, with λ ∝ D R / H , and then collapse at fields distributed around H c ∝ D R 4 / 3 . The concentration of the skyrmions goes down with the field as exp [ − 1 2 ( H / H c ) 3 / 2 ] .


Introduction
Studies of systems with quenched randomness go back to the seminal paper of Imry and Ma (IM) [1] who observed that static disorder, no matter how weak, destroys the long-range order in a system with any continuous-symmetry order parameter in less than four dimensions. The resulting state has received the name of IM domains. The IM argument has been widely used in condensed matter in application to random magnets [2], arrays of magnetic bubbles [3], magnetic flux lattices in superconductors [4], superfluid He-3 [5], superconductor-insulator transition [6], etc. Later analytical studies based upon renormalization group and replica-symmetry breaking methods [7][8][9] suggested that the ordering can be more robust against static randomness, leading to a vortex-free Bragg glass state with an algebraic decay of correlations. Numerical work on systems with quenched randomness [10][11][12], from early on, uncovered strong non-equilibrium effects.
More recently it was found [13] that onset of the irreversible (or glassy) behavior in the presence of a static random field is determined by topology. It depends on the absence or presence of topological defects in the model, as well as on whether the defects are singular or not. The latter is determined by the relation between the dimensionality of space d and the number n of the components of the order parameter, e.g., a fixed-length spin field S r ( ). At n d 1 > + , when topological defects are absent, the behavior is reversible and is characterized by the exponential decay of correlations in accordance with the IM argument. At n d  the presence of topological defects leads to strong metastability and glassy behavior. At n=d+1, when topological defects are nonsingular, the system possesses weak metastability and exhibits a narrow hysteresis loop. This is the case of a twodimensional ferromagnet in which topological defects are skyrmions. Magnetic skyrmions have been subject of intensive recent theoretical and experimental research [14][15][16] due to their promise for developing topologically protected data storage. Mathematically they emerge as solutions [17][18][19]  antiskyrmions, while greater Q | | describe more complex topological objects. In a pure exchange model on a lattice, skyrmions collapse due to violation of the scale invariance by the lattice [20]. Anisotropy, dipole-dipole interaction (DDI), and magnetic field can stabilize significantly large magnetic bubbles with skyrmion topology [21][22][23], while stability of small skyrmions requires other than Heisenberg exchange coupling or a noncentrosymmetric system with large Dzyaloshinskii-Moriya interaction [24][25][26][27][28]. It was shown that quantum fluctuations help stabilize skyrmions as well [29].
In this paper we show that weak quenched randomness provides stabilization of skyrmions in a generic 2D Heisenberg model on a lattice. To make our model closer to reality we consider random anisotropy (RA) rather than random field. In a film geometry the DDI effectively provides an easy-plane anisotropy. To elucidate the effect of the RA we shall assume that the RA energy and Zeeman interaction of spins with the external field applied perpendicular to the film are large compared to the DDI. In this case the effect of the DDI can be neglected. The physical mechanism of the stabilization of skyrmions in such a model is this. Fluctuations of the RA create randomly oscillating energy landscape of the average depth that scales as the square root of the number of spins, N 2 l µ , on a spatial scale λ. Consequently a skyrmion of size λ finds a potential minimum of depth proportional to λ. The Zeeman energy of the skyrmion in the field opposite to its magnetization scales as λ 2 . This provides the minimum of the energy on λ. Discreteness of the crystal lattice adds negative exchange energy proportional to 1/λ 2 [20]. This makes the skyrmions eventually collapse on increasing the field. We begin with the analytical approach and then confirm our picture by numerical studies of spin lattices containing up to 4×10 7 spins. Properties of the skyrmion glass that we observe in the numerical experiment agree well with our analytical findings.

Analytical theory
We consider a constant-length spin field S(r) with the energy where n(r) is a unit vector of the RA and H is the magnetic field applied in the z-direction. Parameters S, α, and β R are related to the parameters of the lattice model with the energy via S=s/a 2 , α=Ja 4 , and D a 2 where a is the lattice spacing. At H=0 the system breaks into IM domains [1] which corresponds to a finite ferromagnetic correlation length [2,19], R a J D f R µ . The argument goes like this. If a significant rotation of the magnetization were to occur on a scale R, the density of the exchange energy would be of order J/R 2 , while fluctuations of the RA energy density, that provide the coherent anisotropy, would be of order D R R a D aR Minimization of the sum of the two energies on R yields the above scaling of R f . This argument should be modified in the presence of topological defects. Such defects in the isotropic 2D exchange model are Belavin-Polyakov (BP) skyrmions. In a continuous spin-field approximation, their energy, due to the scale invariance, does not depend on the skyrmion size λ. Notice that any interaction that breaks the scale invariance, including the RA, destroys the BP solution [30]. However, when a weak interaction is turned on, the skyrmion may still be considered an approximate solution. It will begin to expand or collapse towards the value of λ that corresponds to the energy minimum. This must apply to the RA and Zeeman energies that are typically small compared to the exchange interaction.
The z-component of the spin-field of the skyrmion located at the origin in the background of spins looking up is given by [19] where we assumed λ = L. Writing E Z as HSπR 2 (λ) one can define the magnetic radius of the skyrmion, In the presence of the RA the system size L should be replaced by the length of order R f or Ja/(sH), whichever is smaller. In what follows we will approximate the magnetic size of the skyrmion and its Zeeman energy by where l is a logarithmic factor that we shall treat as a constant.
Statistical fluctuations of the RA result in a coherent magnetic anisotropy on the scale R(λ). As we shall see its strength is proportional to the square root of the number of spins on that scale, which is linear on R(λ). The corresponding pinning energy is given by the statistical dispersion of the anisotropy energy in the region of size r d 2 ò matching the skyrmion size R(λ), To calculate this, it is convenient to switch from the fixed-length anisotropy vector n to a Gaussian distributed vector h r n r µ ( ) ( ) of the average strength h 2 for the average fluctuation of the RA in the area of size r d 2 ò . We use the fact that S(r) in a skyrmion changes on a scale R(λ) that is much greater than the scale, a, of the RA change, which decouples the averaging over h(r) from S(r). Replacing r d 2 ò with π R 2 (λ), we have has been computed in [20]. Introducing it is convenient to write the total energy in the form where H H H c = , c l l l = , and Equation (12) describes three competing effects. The first term is a positive Zeeman energy that tends to shrink the skyrmion because its magnetic moment is opposite to the field. The second term is a negative pinning energy due to the fluctuations of RA. It tends to expand the skyrmion. The third term is a negative lattice energy that tends to collapse the skyrmion. The dependence of the total energy on λ is shown in figure 1. It

Introducing p E E A A
º -á ñ, one can write the normalized distribution for the fluctuating factor p as , applied with probabilities γ and 1 g respectively; γ playing the role of the relaxation constant. High efficiency of this method for glassy systems under the condition 1 g  has been demonstrated by us in the past, see, e.g., [13]. The software used was Wolfram Mathematica that allows compilation (including usage of an external C compiler that doubles the speed) and parallelization. The main operating computer was a 20-core Dell Precision T7610 Workstation. The largest-scale computation, that lasted a few days, has been done on a square lattice containing 6400×6400 spins. Skyrmions have been counted by detecting spots where m z reaches the value 1 and by computing topological charge around such spots.

Numerical results
We begin at H = 0 and apply the field in the negative direction, in small steps, each time minimizing the energy of the system. Topological charge, Q, computed on the lattice by using discretized form of equation (1), is very close to the integer. It is typically small and is due to the statistical difference between the numbers of skyrmions and antiskyrmions. In counting them we do not distinguished between the two. The total number of skyrmions and antiskyrmions is found by looking at the maxima of the z-component of the magnetization, m z . At H=0 skyrmions are comparable in size to the IM domains and are difficult to identify. When the field is applied, the well-separated compact skyrmions with the magnetization opposite to the background emerge, see figure 2. Most of the observed topological defects are simple skyrmions and antiskyrmions with Q=±1, while a small fraction are biskyrmions. At small fields concentration of skyrmions is high. As the absolute value of the field increases, skyrmions shrink and begin to collapse. Correlation between the magnetization maxima and the density of the topological charge, confirming that we are dealing with skyrmions, is illustrated in figure 3.
The average skyrmion size in a weak field is independent of H, see figure 4, in accordance with our expectation that it is determined by R f ∼ (J/D R )a at H=0. As the field continues to increase in the negative z-direction, the skyrmions go down in size and become less abundant. The average skyrmion separation d shown in figure 4 is defined as f 1 S , where f S is concentration of skyrmions, that is, their number per spin. On decreasing the field back to zero no new skyrmions are formed but the existing skyrmions grow in size to adjust to the pinning potential. The sizes of skyrmions in the numerically obtained skyrmion forest were evaluated by computing second derivatives of m z at the maxima and comparing it with the known profile of the BP skyrmion.
The dependence of the average skyrmion size on D H R | | on increasing the field is shown in figure 5. In accordance with analytical results, it is close to linear in the intermediate field range, saturates at small fields, and shows collapse of the smallest skyrmions at high fields. Linear part of the graph allows one to extract the parameter l. In accordance with expectation, it has a weak (apparently logarithmic) dependence on D R . Comparing the slope in figure 5 with equation (14) one obtains l 0.4 » . Note that at very small D/J (large R f ) the finite size of the system begins to affect numerical results.