Robust Formation of Ultrasmall Room-Temperature Neél Skyrmions in Amorphous Ferrimagnets from Atomistic Simulations

Neél skyrmions originate from interfacial Dzyaloshinskii Moriya interaction (DMI). Recent studies have explored using thin-film ferromagnets and ferrimagnets to host Neél skyrmions for spintronic applications. However, it is unclear if ultrasmall (10 nm or less) skyrmions can ever be stabilized at room temperature for practical use in high density parallel racetrack memories. While thicker films can improve stability, DMI decays rapidly away from the interface. As such, spins far away from the interface would experience near-zero DMI, raising question on whether or not unrealistically large DMI is needed to stabilize skyrmions, and whether skyrmions will also collapse away from the interface. To address these questions, we have employed atomistic stochastic Landau-Lifshitz-Gilbert simulations to investigate skyrmions in amorphous ferrimagnetic GdCo. It is revealed that a significant reduction in DMI below that of Pt is sufficient to stabilize ultrasmall skyrmions even in films as thick as 15 nm. Moreover, skyrmions are found to retain a uniform columnar shape across the film thickness due to the long ferrimagnetic exchange length despite the decaying DMI. Our results show that increasing thickness and reducing DMI in GdCo can further reduce the size of skyrmions at room temperature, which is crucial to improve the density and energy efficiency in skyrmion based devices.

www.nature.com/scientificreports www.nature.com/scientificreports/ the compensation temperature 36 . With near zero magnetization and near compensated angular momentum 37 , the skyrmion Hall effect is vastly reduced [18][19][20]23 , and the skyrmion velocity is predicted to be maximum near the compensation point of angular momentum 37 . Recently, all-optical helicity dependent ultrafast switching has been demonstrated in RE-TM alloys using a circularly polarized laser [38][39][40][41][42][43][44][45] . This gives an additional tool to control spins in device structures. Indeed, RE-TM alloys have begun to draw interest in the field of skyrmion research. Large skyrmions of ~150 nm have been observed in Pt/GdFeCo/MgO 23 , and skyrmion bound pairs have been found in Gd/Fe multilayers 24 . Recently, small skyrmions approaching 10 nm were found in Pt/GdCo/TaO x films 25 . Understanding these results will enable further tuning to reduce the size of the skyrmions. To guide experiments, numerical simulation has served as an important tool, especially for complex systems such as RE-TM alloys 40,[46][47][48][49][50] . Several methods, such as atomistic Landau-Lifshitz-Gilbert (LLG) algorithm 40,[46][47][48][49] and micromagnetic Landau-Lifshitz-Bloch (LLB) algorithm 49 have been employed to provide deeper understanding of the magnetic properties in RE-TM alloys.
In this study, an atomistic LLG algorithm 40,[46][47][48] is employed to investigate the properties of skyrmions in GdCo with interfacial DMI. Although the sign of DMI at FM/heavy metal interface is well studied [51][52][53][54][55][56][57] , the sign of DMI involving a FiM remains complex and is rarely discussed. Here, we consider two scenarios for the DMI between Gd and Co (d Gd-Co ). First, DMI between the antiferromagnetic (AFM) pair is set to the same sign as DMI between ferromagnetic pair, i.e. d Gd-Co > 0. Second, the case of d Gd-Co < 0 is considered. The latter appears to be favored by the sign of AFM interaction 27 . Moreover, to incorporate DMI being an interfacial effect, an exponentially decaying DMI is utilized. Simulation results find that with a switched DMI sign, near 10-nm skyrmions remain robust in GdCo films as thick as 15 nm at room temperature. Through numerical tomography maps, we find that skyrmions at room temperature are distributed as a near uniform column in thicker samples, despite a spatially decaying DMI.

Results and Discussion
We will now begin to investigate if ultrasmall skyrmions in GdCo can survive an exponential DMI reduction over thick sample sizes. To incorporate the amorphous nature of GdCo, we employ an amorphous structure of RE 25 TM 75 from ab initio molecular dynamics calculations, as shown in Fig. 1. As shown in Fig. 2, at 300 K, the magnetization of amorphous Gd 25 Co 75 is 5 × 10 4 A/m, and it has a compensation temperature near 250 K. We begin with an exponential DMI decay away from the interface, as shown in Fig. 3. The DMI value discussed herein is the interfacial DMI D 0 . The decay length of DMI is based on both previous simulations and experiments. DMI calculation in Co/Pt interface finds a significant decrease in DMI beyond the second layer of Co from the Pt interface 30 . Experimental results also find similar decay in Co-Alloy/Pt interface [52][53][54] . In amorphous GdCo, we adapted the "second layer" as decay length for direct comparison with these findings.
A range of interfacial DMI values, from d Co-Co = 0.1 × 10 −22 J to d Co-Co = 2.0 × 10 −22 J (D = 0.12 to 2.38 mJ/m 2 ), and three thicknesses, 5 nm, 10 nm, and 15 nm, are considered. Only those show skyrmions are shown herein. To shorten computational time, thicker samples of 10 nm and 15 nm are simulated using a 5 nm thick sample by conserving DMI energy density across the film. To conserve the total DMI energy, a faster decay is employed in a 5 nm sample to simulate 10 nm and 15 nm thick samples to keep the sum of DMI energy to be the same. To check the validity of this simplification, we have compared the results of 10 nm thick samples and 5 nm thick samples with faster DMI decay to verify that the two sets of samples produce identical results. First, we consider www.nature.com/scientificreports www.nature.com/scientificreports/ two scenarios for the sign of d Gd-Co, as both + and − signs have been reported in antiferromagnetically coupled systems 58,59 . Figure 4 shows the color maps of equilibrium spin configurations at 300 K for both d Gd-Co > 0 and d Gd-Co < 0. For the case of d Gd-Gd , d Co-Co > 0 and d Gd-Co > 0, the simulation with d Co-Co = 0.25 × 10 −22 J, d Gd-Gd = 2.96 × 10 −22 J and d Gd-Co = 0.86 × 10 −22 J corresponds to an average DMI of D = 0.21 mJ/m 2 . The value of d Gd-Gd and d Gd-Co is calculated from d Co-Co by multiplying the ratio of Gd moment μ Gd over Co moment μ Co . Further discussion in the supplementary material shows that for a given average DMI, energy minimum and thus skyrmion size is independent of how each DMI term varies. Eq. 3 shows the formula used for converting atomistic DMI to average DMI for Gd x Co 1-x . where n is the average number of nearest neighbors around all atoms, − n A B is the average number of atoms A that are nearest neighbors to atom B, − r A B is average distance between atoms A and nearest neighboring atom B. The π 2 factor comes from averaging of the cross product × s s i j in DMI energy. For 5 nm GdCo, with d Gd-Co < 0, d Co-Co < 0.25 × 10 −22 J, only ferrimagnetic states are observed. At d Co-Co > 1.0 × 10 −22 J, skyrmions are elongated due to boundary effect in the simulation or stripes states are observed. The range of DMI, where skyrmions are found, is smaller compared to calculation by Cort et al. 21 . This is due to a reduction in anisotropy and exchange stiffness in GdCo. With less DMI energy required to create skyrmions, smaller DMI value is needed to create skyrmions and stripes in FiM. Furthermore, experiment results have measured DMI value greater than 1 mJ/m 2 only at ordered FM/heavy metal interface [50][51][52][53][54][55][56] . The DMI value at  www.nature.com/scientificreports www.nature.com/scientificreports/ amorphous FiM/heavy metal remains unknown. Due to disorder nature of amorphous materials, the DMI value in amorphous FiM can be much smaller than the DMI value observed in ordered FM.
As shown in Fig. 4, with ferromagnetic DMI (d Gd-Gd and d Co-Co ) that are positive, two scenarios of AFM DMI (d Gd-Co ) are considered. At 300 K, in all thicknesses, larger DMI is needed to form skyrmions with positive d Gd-Co than with negative d Gd-Co . In 5 nm sample, D = 0.55 mJ/m 2 is needed to stabilize skyrmions with d Gd-Co > 0. In comparison, with d Gd-Co < 0, a smaller DMI of D = 0.21 mJ/m 2 is needed to stabilize skyrmions. Similar behaviors are also found in 10 nm and 15 nm samples. With d Gd-Co > 0, the smallest skyrmions are found at D = 1.26 mJ/m 2 in 10 nm sample and D = 2.31 mJ/m 2 in 15 nm sample. On the other hand, with d Gd-Co < 0, the smallest skyrmions are found at D = 0.84 mJ/m 2 in 10 nm sample and D = 1.68 mJ/m 2 in 15 nm sample.
To understand such intriguing behavior in a FiM, the in-plane spin configurations and the chirality of the skyrmion wall are investigated. Figure 5 summarizes the chirality of the skyrmion walls in the Co sublattice. Using d Gd-Gd , d Co-Co > 0 and d Gd-Co < 0, in the Co sublattice, the spins are turning in counter-clockwise direction across the skyrmion wall. For Gd sublattice, the spins in the skyrmion wall are also turning counter-clockwise. This can be explained by the DMI in the system. AFM couplings between Gd and Co align the spins of Gd and Co in nearly antiparallel directions, except a small canting due to the presence of DMI. With d Gd-Gd and d Co-Co > 0, turning counter-clockwise is energetically favorable. However, with d Gd-Gd , d Co-Co > 0 and d Gd-Co > 0, the chirality of the simulated skyrmion wall is found to be opposite. The DMI torque between the AFM pairs now opposes the DMI torques within each sublattice. In the presence of a stronger inter-sublattice DMI torque, the spins in each sublattice now turn clockwise across the skyrmion wall. www.nature.com/scientificreports www.nature.com/scientificreports/ To better illustrate the change in chirality, the total DMI energies between Co-Co, Gd-Gd and Gd-Co are computed using the equilibrium configurations at 0 K. Table 1 summarizes the sign of the total DMI energies for different nearest neighbor pairs. With d Gd-Gd , d Co-Co > 0 and d Gd-Co > 0, spins are turning counter-clockwise. With this configuration, the total DMI energy between Gd-Gd pair E DMI (Gd-Gd) and Co-Co pair E DMI (Co-Co) are negative, and the total DMI energy between Gd and Co pair E DMI (Gd-Co) is also negative. This means that with d Gd-Gd , d Co-Co > 0, it is energetically favorable for spins to turn counterclockwise. On the other hand, with d Gd-Gd , d Co-Co > 0 and d Gd-Co > 0, spins are revealed to turn clockwise from the simulated configurations. As a result of the sign change in chirality, E DMI (Gd-Gd) and E DMI (Co-Co) become positive. On the other hand, E DMI (Gd-Co) remains negative, because both chirality and d Gd-Co changes sign. This implies that it is energetically favorable for Gd-Co pair to turn clockwise across, but it is energetically unfavorable for Gd-Gd and Co-Co pairs to do so. In other word, AFM DMI d Gd-Co is able to overcome ferromagnetic DMI d Gd-Gd and d Co-Co, resulting in energy favorable configurations for Gd-Co pairs. To summarize, in a FiM, if the DMI of ferromagnetic pair and AFM pair have the same sign, a cancellation of DMI occurs because it is preferable for a ferromagnetic pair to turn in the opposite direction of an AFM pair. No cancellation occurs if the DMI of ferromagnetic pair and AFM pair have the opposite sign. These also explain the differences in size of skyrmion between d Gd-Co < 0 and d Gd-Co > 0. The d Gd-Co < 0 scenario has larger skyrmions because both ferromagnetic and AFM pairs are contributing to the formation of a skyrmion, which means the DMI effect is stronger overall.
To investigate the minimal size of room temperature skyrmions in GdCo, D-K phase diagrams with exponentially decaying DMI at 300 K are simulated for 5, 10 and 15 nm GdCo films. In this section, we focus on the d Gd-Co < 0 scenario. Since energy barrier is a function of exchange stiffness and thickness 18 , the minimal skyrmions size found in d Gd-Co < 0 scenario can also apply to d Gd-Co > 0 scenario, except a larger DMI is required. For each thickness, anisotropy ranges from 0.05 × 10 5 J/m 3 to 4 × 10 5 J/m 3 are investigated. Experimentally, GdCo has anisotropy in the order of 10 4 J/m 3 25,36 . For DMI, larger interfacial DMI is explored in thicker samples, because as thickness increases, the average DMI decreases, and larger interfacial DMI is needed to stabilize skyrmions. In 5 nm samples, interfacial DMI of 0 to 2 mJ/m 2 , which corresponds to d Co-Co of 0 to 2.38 × 10 22 J, are investigated. Figure 6(a) shows the D-K phase diagram of 5 nm GdCo at 300 K. In 5 nm GdCo, skyrmions range from 12 nm to 40 nm are stabilized in the simulated range of interfacial DMI and anisotropy. Lines of 15 to 30 nm indicate the size of skyrmions at various DMI and anisotropy. As DMI decreases or anisotropy increases, skyrmions become smaller and eventually collapse into FiM states. At the opposite side of D-K diagram, with large DMI and small anisotropy, skyrmions larger than 40 nm becomes elongated or collapsed due to the boundary of the simulation space (50.7 nm × 50.7 nm). This elongation of skyrmions was also seen earlier in Fig. 4 at large DMI values. Overall, for a given anisotropy, as interfacial DMI increases from 0 to 2.0 mJ/m 2 , the equilibrium configuration goes from FiM to skyrmions, then to stripes. For a fixed DMI, as anisotropy increases, size of skyrmions decreases, and finally, skyrmions collapse into FiM states. These behavior of skyrmions in FiM GdCo as a function of DMI and anisotropy is the same as what has been observed in a ferromagnet 17,18 .    25,36 , which is in the order of 10 4 J/m 3 . However, the interfacial DMI is less than what has been observed at a Pt interface. Ab-inito calculation has found Interfacial DMI of up to 12 mJ/m 2 is reported at an ideal Pt/Co interface 30 . On the other hand, the interfacial DMI measured in Co/Pt and other Co-alloy/Pt films are around 1.2 to 1.5 mJ/m 2 52-54 . Thus, some reductions of DMI from that of Pt are needed to www.nature.com/scientificreports www.nature.com/scientificreports/ experimentally obtain ultrasmall skyrmion in 5 and 10 nm GdCo films. Reduction of DMI can be obtained by sandwiching GdCo between two Pt layers with one Pt layer being diluted by other elements. Since GdCo is amorphous, we have more flexibility of tuning the underlayer and the capping layer of a multilayer sandwich. With its intrinsic anisotropy and flexibility, GdCo films are promising materials to obtain ultrasmall skyrmions at room temperature through DMI tuning.

Scenario E DMI (Gd-Gd) E DMI (Co-Co) E DMI (Gd-Co)
For device applications, especially in thicker films, we will also need to consider the growth of skyrmions away from the interface. With decaying DMI away from the interface, spins at the top of a thicker sample experience effectively zero DMI. Without DMI, one might expect spins near the top to align parallel for FM neighbors and antiparallel for AFM neighbors, and skyrmions to disappear far away from the interface. If skyrmions collapse far away from the interface, the reliability of such memory devices would be vastly reduced. To investigate whether skyrmions remain robust in thicker samples, a numerical tomography is employed to image simulated ultrasmall skyrmions at 300 K. Figure 7 shows the numerical tomography plot of a ultrasmall skyrmion in 10 nm GdCo. This skyrmion corresponds to D = 0.84 mJ/m 2 and K = 0.3 × 10 5 J/m 3 . The same skyrmion was shown in Fig. 4(b) and as the smallest skyrmions (Star Symbol) at K = 0.3 × 10 5 J/m 3 in Fig. 6(b). In the 3D plots at the center of Fig. 7, colors are made to be somewhat transparent to reveal the skyrmions structure near the center. For Co sublattice, red to orange color shows that most of the spins are pointing down. A region of green and blue that appears near the center corresponds to the simulated skyrmion at 300 K. As evidenced by the columnar distribution of blue color, the skyrmion retains a uniform columnar growth from the bottom to the top. Columnar distribution of skyrmion is also found in Gd sublattice, where a column of red is distributed uniformly from the bottom to the top. This feature can be understood in terms of the large magnetic exchange length >30 nm due to the low magnetization in the ferrimagnet.
To further demonstrate the uniform columnar distribution of skyrmion, in-plane and out-of-plane cross sections of the skyrmion are also plotted in Fig. 7. On the left of Fig. 7, in-plane cross section of spin configuration within 0.5 nm of the interface and 0.5 nm of the top are mapped. The skyrmions at the interface and near the top have identical size and shape. Compare to the mapping of spin configurations in Fig. 4(b), size of the skyrmion remain the same. This shows that the size of skyrmions remain the uniform throughout a sample. On the right side of Fig. 7, out-of-plane cross sections are shown for Gd and Co sublattices. The blue color in Co sublattice and the red color in Gd sublattice correspond to the center of the skyrmion. For both sublattice, out-of-plane cross sections show a columnar distribution of skyrmion from the bottom interface to the top. These results provide important evidences that skyrmion remain robust through a thicker sample, and further support of using thicker GdCo samples to increase skyrmion stability at room temperature. www.nature.com/scientificreports www.nature.com/scientificreports/

Conclusions
Using atomistic stochastic LLG simulations, ultrasmall skyrmions are shown to be stable at room temperature in ferrimagnetic GdCo. Despite the rapid decay of Dzyaloshinskii Moriya interaction (DMI) away from the interface, a realistic range of DMI values is seen to stabilize skyrmions in GdCo films as thick as 15 nm irrespective of the sign of DMI between antiferromagnetic coupled Gd and Co, Furthermore, the low DMI values needed to form ultrasmall skyrmion in GdCo indicate opportunity for designing magnetic materials to host ultrasmall Neel skyrmions. Through tomography of an ultrasmall skyrmion in 10-nm thick GdCo film, it is discovered that the skyrmion assumes a columnar configuration that extends uniformly across the film thickness despite having near zero DMI far away from the interface. These findings argue for using thicker magnetic films to host ultrasmall skyrmions, providing an important strategy for developing high density and high efficiency skyrmion based devices.

Methods
The classical atomistic Hamiltonian H in Eq. (1) is employed to investigate magnetic textures in amorphous FiMs.
where s s , i j are the normalized spins and µ µ , To find the ground state, the spins are evolved under the following stochastic Landau-Lifshitz-Gilbert (LLG) equation, where γ is the gyromagnetic ratio, α is the Gilbert damping constant, H eff is the effective field, ξ is the Gaussian white noise term for thermal fluctuations and M s is the saturation magnetization. The parameters used in our simulation are listed in Table 2. Exchange couplings J ij are calibrated based on Ostler et al. 49 to maintain the same Curie temperature and compensation temperature for a given compensation. At 300 K, the magnetization of Gd 25 Co 75 is 5 × 10 4 A/m. Anisotropy energy is determined based on Hansen et al. 36 To incorporate the amorphous short range order, an amorphous structure of a 1.6 nm × 1.6 nm × 1.6 nm box containing 250 atoms is generated from ab initio molecular dynamics calculations by Sheng et al. 60 . The composition used in the simulation is Gd 25 Co 75 . Figure 1 shows a plot of RE and TM atoms in the amorphous structure. For a 4.8 nm thick sample, replicas of this box (32 × 32 × 3) are placed next to each other to expand the simulated sample to 50.7 nm × 50.7 nm × 4.8 nm and 768000 atoms. On average, we find that each Co atom has 6.8 Co neighbors and 4.1 Gd neighbors, while each Gd atom has 11.7 Co neighbors and 3.5 Gd neighbors. We have also employed a FCC lattice to study skyrmions in GdCo. We found that with the same compensation temperature and magnetization, a larger DMI is needed to stabilize skyrmion in a FCC lattice structures than amorphous structure. This is because the overall effectiveness of DMI is affected by the structure. Only results using the amorphous structure are shown herein.
In the simulations, the initial states are skyrmion of 20 nm based on the 2-pi model 18 . Various initial states, includes random initial states and 10-30 nm skyrmions, have been tested and found to produce the same final states. The size of skyrmions are defined as the diameter for which M z = 0. Since skyrmions are not perfectly symmetric, size of skyrmion is the average diameter.

Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.