Mapping magnetization states in ultrathin films with Dzyaloshinskii–Moriya interaction

To succeed in the next generation of magnetic storages, nanomagnetic structures must be engineered to allow modulation of magnetization amplitude and spatial period. We demonstrate computationally how magnetic nanostructures states (domains with narrow wall, skyrmions, spin spirals, conical spin spirals, and in-plane magnetization configuration) can be designed in ultrathin films with a Dzyaloshinskii–Moriya interaction (DMI) by adjusting two material parameters: perpendicular magnetic anisotropy characterized by the quality factor Q, and reduced DMI constant Δ . For a broad range of Q and Δ parameters, the magnetization states are mapped in (Q, Δ ) diagrams and characterized by the periodicity (p) of spatial distribution of magnetization and the mean value of the square of an out-of-plane normalized magnetization component 〈 m z 2 〉 . We explain the transitions between different magnetization states and describe phase transition curves in terms of vanishing of the domain wall energy for both 0 ≤ Q < 1 and Q > 1 cases. We show that by changing Δ and/or Q and approaching the phase transition curves, the discontinuous transitions accompanied by jumps of both period and amplitude take place.


Abstract
To succeed in the next generation of magnetic storages, nanomagnetic structures must be engineered to allow modulation of magnetization amplitude and spatial period. We demonstrate computationally how magnetic nanostructures states (domains with narrow wall, skyrmions, spin spirals, conical spin spirals, and in-plane magnetization configuration) can be designed in ultrathin films with a Dzyaloshinskii-Moriya interaction (DMI) by adjusting two material parameters: perpendicular magnetic anisotropy characterized by the quality factor Q, and reduced DMI constant D. For a broad range of Q and D parameters, the magnetization states are mapped in (Q, D) diagrams and characterized by the periodicity (p) of spatial distribution of magnetization and the mean value of the square of an out-of-plane normalized magnetization component á ñ m z 2 . We explain the transitions between different magnetization states and describe phase transition curves in terms of vanishing of the domain wall energy for both  < Q 0 1 and Q>1 cases. We show that by changing D and/or Q and approaching the phase transition curves, the discontinuous transitions accompanied by jumps of both period and amplitude take place.
Magnetism on nanometer scale is a subject of great interest in view of both the development of fundamental knowledge and its possible applications for next generation high density magnetic recording. Magnetic skyrmions were derived first in 1989 [1] and they were being intensively investigated [2]. They attract attention due to their huge application potential as a carrier of information in racetrack memory [3,4]. Among others, spectacular successes have been obtained, as controlled creation, shifting and annihilation of skyrmions in multilayer metallic systems at room temperature, using for manipulation relatively low magnetic fields and electrical currents. Large efforts have been devoted to engineer such metallic multilayer materials, in which the Dzyaloshinskii-Moriya interaction (DMI) exists as large as possible, being necessary to stabilize skyrmions at room temperature [5,6]. Besides skyrmions, another interesting magnetic structure of the same class-enabled by DMI, is spin spiral, with one-dimensional (1D) chiral modulations of spins [7,8].
Strength of DMI is defined by the constant D from E DMI =D·(S 1 ×S 2 ), where E DMI is the energy of DMI, D is the Dzyaloshinskii-Moriya vector, and S 1 and S 2 are two coupled spins. Changes of D is possible by proper adjustment of magnetic layer interfaces [4,9]. A key role of out-of-plane magnetic anisotropy, K, in establishing magnetization spatial phases has been described [10,11]. Different methods to change the magnetic anisotropy are known: (i) choosing a multilayer structure during sample deposition (ii) sample post-growth modification of either by light pulses [12] or by ion irradiation [11,13], (iii) application of the electric field [14][15][16]. Geometrical constraints of ultrathin layers can also drive the stability and evolution of magnetic states [17].
A phenomenological quasi-classical continuum approximation can be used to calculate the skyrmion energies and length scales [7]. In particular, it was shown that the spontaneous skyrmion lattice ground states may exist in thin films, where a lack of space inversion symmetry leads to chiral interactions. Recently, the stability of small radius skyrmion in an infinite film as a function of D and perpendicular magnetic field was calculated [18]. During recent years several aspects of DMI were studied experimentally, like an influence of magnetic anisotropy on skyrmions, in continuous ultrathin films [19] and cylindrical nanodots [20], Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. temperature evolution of spin spirals [8], magnetic domain structures driven by the interfacial DMI in perpendicularly magnetized exchange biased systems [21], or by the DMI modulated by electric field [22]. A list of materials that are known to host skyrmion lattices or single skyrmions and zoo (topological spin configurations) of skyrmion's types were recently systematized [23].
Since both K and D parameters can be easily driven, e.g. by the thicknesses and/or chemical composition of surrounding layers during deposition, as well as by other post-growth manipulations, the main goal of the present paper is to perform micromagnetic simulations and systematic mapping of magnetization states in ultrathin magnetic layers for a wide range of both anisotropy and DMI magnitude, and to study their evolution with an externally applied magnetic field. We use the following reduced material parameters throughout, unless specified otherwise. Magnetic anisotropy is expressed by the quality factor Q, i.e. the ratio of the magnetic anisotropy and demagnetization energy: The DMI constant is also used in a reduced form where D 0 is defined as follows: Micromagnetic simulations have been performed with the OOMMF software [24]. The total system energy = consisted of the exchange interaction, uniaxial anisotropy energy, Zeeman energy, demagnetization (stray field) energy, and DMI, realized by the extension package [25]. In the limiting case without DMI, the energy (equation (5)) describes the classical stripe domains, isolated bubble domain, bubble lattices, and sinusoidal domain structures [10,26]. On the other hand, including DMI and neglecting the demagnetization energy, leads to appearance of chiral modulated states [2,27]. Similarly to our previous simulations for cobalt [10,11], we took the constant values of the saturation magnetization M S =1420 kA m −1 and the exchange constant A=1.3×10 −11 J m −1 . Thus, the value of D 0 can be calculated with equation (3), D 0 =5.16 -mJ m 2 . The thickness of magnetic film was also kept constant, equal to 1.86 nm (which is a spin reorientation transition (SRT)thickness in an infinite Co film, according to our previous simulations [10]). As the independent parameters of the problem, the quality factor Q was changed between 0 and 3, and reduced DMI constant D was changed between 0 and 2.
The simulations were carried out in either 1D or two-dimensional (2D) mode. A choice of 1D/2D simulation mode was made depending on the studied objects. The 1D mode was found to be satisfactory for all types of stripe-like structures. Although the domains appeared then as un-natural infinite straight lines, the crucial characterizing parameters were found to be consistent with the ones obtained from 2D mode, as was tested for a series of Q and D combinations. An undoubtful advantage of 1D mode was that bigger systems could be modeled, regarding the limitations of memory and computation time. The 2D mode was used, when a planar distribution of magnetization was important-particularly for skyrmions and their transformations. In the 1D mode, the simulated sample consisted of a chain of cuboid cells, each of 2×3000×1.86 nm 3 , distributed along the x-axis and forming a cuboid structure of either 3000×3000×1.86 nm 3 or 10000×3000×1.86 nm 3 in size, depending on the expected period of simulated stripes. The larger sample has been used for simulation in case of large Q and small D parameters. In the both cases, the sizes of the simulation cell in the xz plane was smaller than the exchange length of cobalt, equal to 3.2 nm. In the 2D mode, the simulated sample consisted of a square array of cuboid cells, each of 2×2×1.86 nm 3 , forming a cuboid structure either 200×200×1.86 nm 3 or 1000×1000×1.86 nm 3 in size, depending on the expected sizes of magnetic objects. The periodic boundary conditions have been switched on along both x and y directions, avoiding boundary effects. The damping coefficient has been set to 0.2. The magnetic field was applied perpendicularly to the sample plane (along z direction). The calculations were terminated when the maximal absolute magnitude of the time derivative of magnetization across all spins dropped below one hundredth degree per nanosecond in the case of 1D mode and half a degree per nanosecond in the case of 2D mode.
Let us discuss the results of micromagnetic simulations performed without an external magnetic field. The results of the 1D mode simulations for the selected values of Q and D, are shown in figures 1-3. As was mentioned, a consistency of these 1D results with 2D mode was controlled for many (a few hundreds) sets of Q and D values, throughout the studied ranges of the both parameters. It is well known that in the absence of the DMI (i.e. D=0), the SRT from homogeneous out-of-plane mono-domain state to homogeneous in-plane magnetization configuration occurs at Q=1. Considering film splitting into a stripe magnetic domain structure, we found a drastic decrease of the domain period, while decreasing Q, until critical value * p =8πl ex 2 /d+2πd (where d is the film thickness) [10], achieved at the critical * Q smaller than 1. Taking this into account, in our current modeling of magnetization distributions for different D, we have selected four representative values of Q far from the SRT and in the SRT vicinity: (i) Q 1 =0, zero uniaxial magnetic anisotropy, (ii) Q 2 =0.7-low anisotropy, (iii) Q 3 =0.98 which is slightly lower than the critical * Q -close to the SRT, and (iv) Q 4 =1.1, relatively high anisotropy.  The bottom row of images in figure 1(a), obtained for D=0 and the selected values of Q, depicts the classical SRT scenario [10]. While decreasing Q, the magnetization distribution evolves: (i) from out-of-plane oriented, large period domain structure with narrow Bloch domain walls (Q 4 =1.1, 4th column, structure B in figure 1(b)) (ii) through the sinusoidal domain structure (magnetization oscillation in the yz plane, structure D in figure 1(b)) with a small critical period close to * p (Q 3 = 0.98, 3rd column) (iii) to the in-plane magnetization state in the case of (Q 2 = 0.7 and Q 1 = 0, the 2nd and 1st columns, structure F in figure 1(b)).
Looking on the panels in columns shown in figure 1(a), from bottom to top, one can see how the D-driven evolution of spatial magnetization distribution undergoes when changing the values of Q.
In case of Q 4 =1.1, while increasing D, one can find the two following steps of magnetization distribution evolution (see the 4th column in figure 1(a)): (i) the transformation of Bloch-type domain walls into Néel-type walls and a decrease of period (for D=0.2, structure A in figure 1(b)); (ii) a further decrease of the domain period and magnetization distribution transformation towards the spin spiral domain structure (for D=0.4, structure C in figure 1(b)). The decrease of the domain period connected with the increase of D value is clearly visible in figure 2(a), where the dependence of the period on D is plotted. The inset to figure 2(b) shows D-driven transformation of m z (x/p), a profile of the out-of-plane magnetization component across the normalized period, from a domain structure with a narrow domain wall (low D) into a spin spiral-like structure (large D). We characterize the simulated magnetization states by á ñ m z 2 , averaged over the whole sample, which is a convenient measure of 'magnetization perpendicularity', e.g. for: the homogenous perpendicular magnetization state á ñ m z 2 =1, out-of-plane domains with narrow walls á ñ m z 2 <1, spin spiral structure á ñ m z 2 ≈1/2, sinusoidal domain structure 0<á ñ m z 2 1/2, and in-plane state á ñ m z 2 =0. The changes of the type of domain structure and corresponding changes of á ñ m z 2 are illustrated in figure 2(b). In case of D=0 and Q 4 =1.1, the domain structure has an m z amplitude component close to unity, where a large period and narrow Bloch walls exist, making á ñ m z 2 ;1. While D increases, the á ñ m z 2 parameter decreases continuously down to 1/2, which is typical for a zero-field spin spiral domain structure (if m z (x) followed the cosine function, á ñ m z 2 would be exactly 1/2). Typical approaches to analyze the magnetization distribution require assumptions about domain walls type and parameters [26], and initially it should be decided what is an applicable model to the given structure. Here, without any preliminary knowledge, the structures can be characterized simply by á ñ m z 2 allowing to qualify one of several possible magnetization states.  figure 1, the bottom panel, structure D in figure 1(b)) and oscillating with x. Sinusoidal domain structure exhibits Bloch-like domain walls, which width is comparable in size with the domain period. While increasing D, one can find the two following steps of magnetization distribution  figure 2(b) it is visible, that the parameter á ñ m z 2 starts from the value slightly above 0 (D=0) for the sinusoidal domain structure, next grows rapidly to values greater than 0.7, reflecting an increase of the out-of-plane magnetization component and finally decreases towards 0.5, what is connected with the decrease of the domain period and transformation of magnetization distribution towards the spin spiral state.
In case of Q 2 =0.7 and D lower than 0.61 the in-plane magnetization state exists-see the 2nd column in figure 1. While increasing D, the two following evolution steps of spatial magnetization distribution undergo: (i) at D=0.61 stripe domains appear with conical spin spiral magnetization distribution (structure E in figure 1(b)), (ii) for larger D, a decrease of domain period and magnetization distribution transformation towards the spin spiral domain structure is observed (see the panel for D=1.2).
A similar D-induced magnetization distribution evolution undergoes in the case of Q 1 =0, see figure 1the first column, but now the stripe domains with conical spin spiral magnetization distribution appear at higher D, namely D=1.07.
For these both last cases of low anisotropy (Q=0 and 0.7), D-induced magnetization distribution evolution is presented in figure 2, too. While increasing D, the domain period decrease is visible in figure 2(a). In figure 2(b) one can see the growing dependence of the á ñ m z 2 parameter on D, starting from zero, which corresponds to the in-plane oriented homogenous state of magnetization, towards a value close to 1/2, which corresponds to the spin spiral state.
Summarizing this part, all images in top row of figure 1(a) present spin spirals, as theoretically predicted for the regime, where the demagnetization energy is negligibly small in comparison with the DMI energy [2,27].
For a zero field and zero effective anisotropy the spin spiral period is given as [27]: The p ss (Δ) dependence is drown as the black thick line in figure 2(a). This line is close to the simulated periods of spin configurations in a wide range of D, except of D close to zero, where the magnetic anisotropy and demagnetization contributions mainly determine the period value. While decreasing D to zero the theoretical spin spiral period D ( ) p ss goes to infinity, the periods of the simulated structures for both Q 3 =0.98 and Q 4 =1.1 go to finite values. This deviation is a consequence of neglecting the demagnetization contribution, while obtaining equation (6). Remarkable is that the all simulated values of the domain period are comparable with the spin spiral period given by equation (6).
An extension of the results presented above is given in figure 3 in the form of D ( ) Q, diagrams of period p ( figure 3(a)) and á ñ m z 2 parameter ( figure 3(b)). To analyze the diagrams let us consider energies of the isolated 180 • domain wall for Q > 1 and Q < 1 where magnetization in domains is directed perpendicularly and inplane, respectively. Three contributions into wall energy should be considered: magnetic anisotropy, demagnetization and Dzyaloshinskii-Moriya.
After the work [28] for out-of-plane domains and Q > 1, the domain wall energy can be ascribed as For in-plane domains and Q < 1, the Bloch wall energy could be calculated in a similar way When the domain wall energy starts to be negative at a sufficiently high value of D (or D), the domains would proliferate in the sample, allowing a domain period to decrease ( figure 3(a)). Equalizing the wall energies to zero in equations (7) and (8), one straight forwardly obtains the equations for the critical lines: For D > D^| | crit , the negative wall energy manifests in the instability of the homogeneous states with respect to chiral modulations [29].
The boundary lines D^( ) || Q crit , are plotted in figure 3 with red lines. Reaching the boundary lines, depending on Q, the two scenarios of the transitions between different magnetization states can be distinguished (see figures 2 and 3) while increasing D: (i) the transition between the in-plane magnetization and spin spiral states and (ii) the transition between the domains with narrow walls and spin spiral domain states. As one can see in figure 3(b), the theoretical boundary curves D ( ) || Q crit and Δ crit⊥ (Q) fit very well with the contours for á ñ m z 2 =0.05 and á ñ m z 2 =0.8, respectively. We will now discuss the results of simulations of magnetization distributions in a presence of a perpendicular magnetic field. Generally, the magnetic state of the sample strongly depends also on magnetic history and initial magnetization distribution. Different zero-field metastable states are possible: stripe domains (spin spirals), as well as lattices of bubbles (skyrmions) and a mixture of these states. Simulations of 2D sample 200×200 nm 2 were performed for Q 4 =1.1, over a wide range of D, enabling analysis of field-driven evolution for different magnetization distributions, from magnetic domains with narrow walls (at D=0), to spin spiral states (D  0). Different scenarios of magnetic history were realized.
First, as an initial magnetization distribution we chose a single skyrmion of small radius (of a few simulation cells) placed at the center of the homogeneously magnetized sample in the presence of a high enough magnetic field. Then the field was decreasing gradually to zero. In this part of the simulations, we took care about skyrmions to be considered as isolated. Due to the use of periodic boundary conditions, a skyrmion can interact with its periodic copies. While a skyrmion is mainly driven by the short-range interactions, the skyrmionskyrmion interaction represents a long-range magnetic dipole-dipole interaction, which decays on the length scale of the order the skyrmion radius. So, if the sample is big enough, the skyrmion radius is at least one order of magnitude smaller than the sample size, and one can neglect an interaction between a skyrmion and its periodic copies.
Decreasing the field down to a certain critical value, the skyrmion size increased up to a certain critical radius. The field-dependence of the skyrmion radius, r sk (B) (for D=1.4 and Q 4 = 1.1) is plotted in figure 4 with open circles. Below a certain critical field the skyrmion is converted rapidly into stripe magnetic domains. That field can be called 'elliptic', B ell , by an analogy to a situation of cylindrical (bubble) domains, whose shapes under decreasing magnetic field changes to elliptic, before regular stripes are formed. In figure 4 the characteristic size of the stripe domains is also plotted with solid circles, as a quarter of period, p 4, to be able to directly compare the sizes of skyrmions and the stripes. Such an approach is explained in the top part of figure 4. Both quantities, p/4 and r sk , describe the distance, at which the magnetization rotates from out-of-plane to inplane orientation. Knowing p(B) and r sk (B), the normalized out-of-plane component m z (B) was calculated, see the blue curve in figure 4. Another critical field and a boundary, which can be revealed by the simulations, is the equilibrium field, B eq at which the energy of the sample with a single isolated skyrmion is equal to the energy of the homogeneously out-of-plane magnetized sample. B eq field is distinguished in figure 4 by the thin vertical line. The quantity of B eq is particularly important, because below that value the energy of a skyrmion is lower than the energy of the homogeneously magnetized sample, thus the skyrmion state is somehow preferred.
The values of critical field B ell , for Q 4 =1.1, and across the studied values of D, are plotted in figure 5 with open blue squares. At field B ell , i.e. just before a conversion from the skyrmion state to the stripe domains state, a skyrmion reaches a maximal radius, r ell , and this quantity was also found to be strongly dependent on D, as presented in figure 5 with red circles. On decreasing D, the radius r ell increases rapidly. That dependence, as well as D ( ) B ell , is cut below D=0.3, because below that value a very large simulated sample and highly timeconsuming simulations would be required. Figure 5 also presents the D-dependence of B eq , using solid squares. As seen from this figure, B eq rapidly increases with D. Such behavior is plausible, since increasing DMI a larger field is required to prevent the system from developing inhomogeneities. Interestingly, B eq has a non-zero value at D=0, however that point should be assigned to a bubble with a Bloch wall rather than a Néel skyrmion.
In continuous models of skyrmions, with the superimposed boundaries conditions for the orientation of magnetization vector versus the distance from the skyrmion's center, θ(0)=π and q ¥ ( )=0, a skyrmion exists at arbitrary high fields. In most experiments skyrmions are stabilized by DMI and an external magnetic field. At increasing fields, the skyrmion size gradually decreases and asymptotically approaches zero [29]. Here, in a discrete model, it is possible to observe a skyrmion collapse. However, this happens when the skyrmion size becomes comparable with the size of a single micromagnetic cell. Further, the collapse field is found to be dependent on the cell size, what seems to make that micromagnetic approach to the study of the collapse field inaccurate. A more correct way should involve spins on an atomic lattice, finding a minimum energy path between two stable configurations [30]. Such approach is beyond the scope of the present work, but here B eq can be considered as a kind of estimation of the collapse field, as it describes a situation when two stable configurations have equal energies.
In summary, the magnetization distributions in ultrathin magnetic films with DMI were simulated in a wide range of magnetic anisotropy and DMI constant. This allowed to map the equilibrium and metastable magnetization states on a D ( ) Q, diagram (figures 2 and 3). We showed that the interplay between uniaxial anisotropy and DMI energy contributions gives rise to a big variety of domain-like magnetization states and isolated skyrmions. The periods of domain-like magnetization distributions (stripe domains, sinusoidal  domains, domains with the Bloch and Néel domain walls, and spin spirals), calculated as functions of the quality factor and the DMI constant, have the same order of magnitude. When changing Q and D, the transitions between either the in-plane state and spin spiral structure or between the in-plane and out-of-plane domain states with narrow domain walls take place. These transitions are accompanied by abrupt changes of the domain period, as well as an out-of-plane magnetization component ( figure 3). To describe these transitions we give simple analytical formulas allowing to define the phase transition boundaries on a D ( ) Q, phase diagram. We also showed that decreasing the perpendicular magnetic field, an isolated skyrmion exhibits the elliptical instability at the magnetic field B ell when the skyrmion radius reaches its maximal value, r ell ( figure 4). While decreasing D, a drastic decrease of both B ell and B eq fields, as well as an increase of the r ell radius were found. The magnetization state in ultrathin film could be chosen (changing Q, D) by proper sample designing (by deposition condition or post-growth treatment) or by reversible changes of Q, D, e.g. by external electric field or temperature. Taken together, our results give a roadmap for engineering stable and metastable magnetic nanostructures states with properties needed for applications in the next generation of magnetic storages.