Micromagnetic Design of Skyrmionic Materials and Chiral Magnetic Configurations in Patterned Nanostructures for Neuromorphic and Qubit Applications

Our study addresses the problematics of magnetic skyrmions, nanometer-size vortex-like swirling topological defects, broadly studied today for applications in classic, neuromorphic and quantum information technologies. We tackle some challenging issues of material properties versus skyrmion stability and manipulation within a multiple-scale modeling framework, involving complementary ab-initio and micromagnetic frameworks. Ab-initio calculations provide insight into the anatomy of the magnetic anisotropy, the Dzyaloshinskii–Moriya asymmetric exchange interaction (DMI) and their response to a gating electric field. Various multi-layered heterostructures were specially designed to provide electric field tunable perpendicular magnetization and sizeable DMI, which are required for skyrmion occurrence. Landau–Lifshitz–Gilbert micromagnetic calculations in nanometric disks allowed the extraction of material parameter phase diagrams in which magnetic textures were classified according to their topological charge. We identified suitable ranges of magnetic anisotropy, DMI and saturation magnetization for stabilizing skyrmionic ground states or writing/manipulating them using either a spin-transfer torque of a perpendicular current or the electric field. From analyzing the different contributions to the total magnetic free energy, we point out some critical properties influencing the skyrmions’ stability. Finally, we discuss some experimental issues related to the choice of materials or the design of novel magnetic materials compatible with skyrmionic applications.


Introduction
One of the most recent magnetic paradigms in information and communication technologies (ICT) is based on the use of magnetic skyrmions. These vortex-like swirling topological defects in the magnetization texture are particularly interesting due to their extremely small size, and low level of energy consumption when manipulated to change their state that can encode either classic, neuromorphic or quantum information. The request for low energy consumption represents a major challenge of next-generation ICT devices, having in view the exponential increase in energy consumption in this area. In 2018, ICT accounted for between 5% and 9% of total electricity consumption worldwide [1]. Widely cited forecast sources suggest an acceleration of the total energy demand in the ICT area in relation to new blockchain applications (e.g., cryptocurrencies), and to emerging neuromorphic and quantum technologies. Large data centers and wireless and wired networks demand more and more energy compared to consumer devices such as televisions, computers, mobile phones. An optimistic prediction expects that in 2030, more than 20% of the global electricity consumption, about 9000 Terawatt (TWh) hour, will come from the ICT sector [1]. This trend is clearly distinguishable when looking at the data provided by Google on the web, regarding the energy consumption for the search query: an average Google search event requests about 0.0003 kWh of energy, the equivalent of turning a 60 W light bulb for about 17 s, 30 searches cost the energy of boiling one litter of water, the total Google search consumption in 2015 was about 5.7 TWh, similar to the total annual energy consumption of San Francisco (5.8 TWh). This trend tripled between 2015 and 2020, the analysis illustrating a clear exponential increase. The increase will be further amplified by the increase in the total number of bits produced and manipulated annually on Earth [2], today about 10 21 . In the next 300 years, the power requested to sustain the digital production is expected to exceed 18.5 × 10 15 W, i.e., the total planetary power consumption, today. Therefore, aside from the existing global challenges related to the environment, climate, population, health, energy, food and security, it is obvious that the energetic efficiency of any device used in ICT represents a major request. New concepts in data storage and the manipulation area must emerge. Spintronics should subscribe to this paradigm and the magnetic skyrmions are particularly suitable for that. They have the advantage of being particularly small, can be manipulated with extremely low levels of energy when storing classic, neuromorphic or quantum information, and ultimately, bring additional functionality driven from magnetic topology issues into consumer-friendly, low-energy nanoscale electronics. Within this context, there is a strong motivation for the theoretical modeling of skyrmionic materials and devices. Our objective is to identify the critical magnetic parameters of materials for hosting and manipulating skyrmions.
In this paper, using a multiscale modeling approach, we address the problematics of skyrmion stabilization versus critical magnetic properties of nanomaterials. Using abinitio calculations on supercells describing realistic/experimental multilayer architectures, we extract the mechanisms and magnitude of some main magnetic properties: magnetic moments and anisotropy (PMA), the Rashba contribution to the interface asymmetric exchange constant DMI, and their response to an external gating electric field. Then, micromagnetic simulations based on Landau-Lifshitz-Gilbert dynamics, led to various magnetic binary phase diagrams of magnetic nanostructures in which the final relaxed state is classified via the topological charge. The micromagnetic simulation results illustrate critical conditions in which skyrmionic states can be written by pulse currents in nano-disks via STT effects. This allows us to identify the DMI-PMA range in which skyrmionic ground states can be stabilized in nano-disks patterned from magnetic multilayers with various sizes and architectures. Furthermore, the results point out some critical parameters which determine the skyrmionic properties (size, stability in temperature, etc.) and illustrate how the skyrmions can be manipulated by external electric field gating or spin-polarized currents. Finally, we discuss some experimental issues related to the choice of materials or the design of novel magnetic materials compatible with skyrmionic applications in classic, neuromorphic and quantum information technologies.

General Issues about Skyrmions
Magnetic skyrmions were expected to be short-time excitations quickly collapsing on linear singularities. However, the seminal work of Bogdanov and Yablonski (1989) demonstrates that they can exist as a metastable configuration with a sizeable lifetime in some low-symmetry condensed matter systems with broken mirroring symmetry [3]. The mechanism responsible for their stabilization is the asymmetric exchange interaction for either bidimensional (2D) or tridimensional (3D) spin localized states. This is the so-called Dzyaloshinskii-Moriya interaction [4] which typically occurs in non-centrosymmetric magnetic materials (bulk DMI), at the interface between a heavy metal and a ferromagnetic thin film (interface DMI) and in some electric field gated ferromagnetic film surfaces (Rashba DMI). Even if at the beginning, the DMI has been only regarded as a higher order effect occurring between ions already coupled by super-exchange (about 10 −2 times smaller than the direct exchange [5], the DMI was further found to be particularly important in ferromagnetic materials with a large spin-orbit coupling (SOC) H SO (2) . This term tilts the magnetic moments away from the collinearity and its magnitude proportional to the SO interaction included in the SO coupling parameter ζ(r). The low dimensionality of magnetic nanostructures brings additional sizeable contributions to the DMI via the associated additional potential gradients and corresponding intrinsic electric fields. At the interface between magnetic and nonmagnetic materials (heavy metals, insulators, etc.) the lack of inversion symmetry leads to a gradient of the lattice potential. Moreover, in multilayered stacks, additive interface effects would provide enhanced DMI and other interesting features. For example, it is possible to tune the sign of the DMI. Because the vector → D ij is parallel to the film plane in magnetic multilayers with perpendicular magnetization, its orientation is defined by the chemical nature and electronic structure effects at the interface between the nonmagnetic and the ultrathin ferromagnet. This leads to an enhancement of the effective magnitude → D ij with decreasing ferromagnetic layer thickness.
Therefore, the magnetic multilayer stacks provide a large degree of freedom in skyrmion stabilization, mainly based on the sizeable interface DMI and the perpendicular magnetic anisotropy (PMA) in ultrathin magnetic films. Within this area, during the past few years, many experimental studies have been performed. They illustrated the successful stabilization of skyrmions from below to above room temperature (RT). The skyrmion stability and its size (core diameter) are controlled by the interplay between the magnitude of the DMI and the dipolar interactions. Most of the studied skyrmionic systems were based on Co, Fe and CoFeB ultrathin ferromagnetic layers, interfaced with Pt, Ta, Ir or W heavy metals and/or MgO insulator [6][7][8][9][10][11][12][13][14][15][16] (see Table 1). Additive effects, by multiple repetitions of the skyrmionic active magnetic stacks (e.g., [Ir/Co(0.6 nm)/Pt] n = 10 ) were found to enhance the thermal stability of skyrmions, even above the room temperature [9]. This represents an extremely important issue for the technological implementation in RT working ICT devices. The skyrmions can be efficiently displaced either by spin-orbit torque [6] or by spin-transfer torque (STT) effects in skyrmion racetrack memories [16] and can be created and annihilated by electric fields [17].
The expected application area of skyrmions is quite broad, from classic spintronic devices such as racetrack memories [16], to neuromorphic spintronics (artificial skyrmionics neurons/synapses) [18] and, as more recently claimed, to quantum spintronics applications based on skyrmionic helicity qubits [19]. A skyrmion racetrack memory is very similar to a domain wall memory [20]. Here, skyrmions are moved by STT of spin-polarized currents instead of domain walls. The presence or absence of a skyrmion, read by a fixed magnetic tunnel junction head placed near the track represents the classic 1 and 0 states. The smaller size of skyrmions compared to domain walls enables larger areal storage densities. In a neuromorphic device, e.g., skyrmionics memristor, the particle-like behavior of skyrmions and their thermal Brownian motion have strong analogies with the neurotransmitter diffusion. Skyrmions accumulating in a "reservoir" could be the spintronic analogues of leaky integrate-and-fire neurons. Moreover, nonlinear resistance changes in magnetic skyrmion systems, originating from the interplay of magnetoresistance and spin or spin-orbit torques on the skyrmions that either move or distort them, can be exploited for unconventional computing [18]. Finally, the skyrmionic quantum bit [19] takes advantage of the quantum information that a skyrmion can store in its helicity: a quantum state being an arbitrary superposition of skyrmion configurations with distinct helicities, uniformly distributed over the Bloch sphere. Moreover, antiferromagnetically coupled Nanomaterials 2022, 12, 4411 4 of 33 skyrmions in synthetic antiferromagnetic multilayers have been recently proposed as platforms for a qubit coupling scheme where the logical states can be adjusted by electric fields [19]. Table 1. Selection of thin film multilayered materials illustrating Néel skyrmion stabilization. We indicate the material (multilayer system), the measured diameter of the skyrmion core, the magnitude of the DMI |D| mJ m 2 , the temperature of the skyrmion stability and the reference of the paper containing the study.

Multilayer System
Diameter of Skyrmion Core (nm) |D| mJ m 2
As magnetic textures, skyrmions resemble vortices but with the z component of the magnetization unit vector m z varying from +1 to −1. In a magnetic nanostructure, the ansatz magnetization profile of a vortex in spherical coordinates, m = (sin θ cos φ, sin θ sin φ, cos θ); r = (r, ϕ) is θ = θ(r), φ(r) = ϕ ± π 2 , and minimizes the volume and edge surface charges: −∇·m = 0 and m·r = 0. The sign ± in φ(r) determines the chirality of the vortex. In a skyrmion the magnetic configuration results from the competition between the exchange, anisotropy and the chiral DMI interaction and can be roughly described by an inexact but useful ansatz profile [22] θ DW (r) = 2arctan sinh(R/w) sinhr/w as a single 360 • domain wall, with some characteristic parameters R, w being the skyrmion size and skyrmion domain wall width, respectively. To further understand the differences between skyrmions in various materials it is important to emphasize the different types of skyrmions and the microscopic origin of their stability. First, the skyrmions can be distinguished by their respective domain walls (DWs), as right-and left-handed Néel and Bloch spin structures, as illustrated by our micromagnetic simulations ( Figure 1).
The main parameters of the skyrmions are: • The topological charge or the winding number Q = |1| that describes how many times the magnetization vector m can be mapped onto a sphere: In saturated magnetic states Q = 0, Q = |0.5|for magnetic vortices, and more complex structures that wrapping the Bloch sphere multiple times Q >|1|.
There is a direct correlation between the winding number and the stability of the spin structures. Theoretically, the Skyrmionic winding number is an integer (positive or negative). The physical manifestation of this topological quantization is that the magnetic skyrmions are topologically protected from decaying to a uniform magnetization that has a winding number equal to zero. As we will show in the sections dedicated to results and discussion, the spins from the core of the skyrmion are protected by an energy barrier driven by the DMI from the reversed spins from the edge.
the helicity. This allows the illustration of the skyrmion's topological analogy with a sphere, the projection on the Bloch sphere of the magnetization vector field being particularly useful for understanding the helicity qubit ( Figure 2).  The skyrmion parameters are linked to the generalized definition of a nanoscale skyrmion magnetization field: ; Φ( ) = + ; with being the helicity. This allows the illustration of the skyrmion's topological analogy with a sphere, the projection on the Bloch sphere of the magnetization vector field being particularly useful for understanding the helicity qubit ( Figure 2). and |0⟩: = −1: | ⟩ = |0⟩+ |1⟩. The representation is issued following our micromagnetic simulations, see next paragraph.
There is a direct correlation between the winding number and the stability of the spin structures. Theoretically, the Skyrmionic winding number is an integer (positive or negative). The physical manifestation of this topological quantization is that the magnetic skyrmions are topologically protected from decaying to a uniform magnetization that has a winding number equal to zero. As we will show in the sections dedicated to results and discussion, the spins from the core of the skyrmion are protected by an energy barrier driven by the DMI from the reversed spins from the edge. Left: Hedgehog representation of a Néel skyrmion magnetization vector field projected on the Bloch sphere. We indicated the vorticity ω, the topological charge Q and the core polarization P. Right: Bloch sphere and the helicity qubit: encoding the quantum superposition of |0 : m z = +1 and |0 : m z = −1: |ψ = α|0 + β|1 . The representation is issued following our micromagnetic simulations, see next paragraph.
There is a direct correlation between the winding number and the stability of the spin structures. Theoretically, the Skyrmionic winding number is an integer (positive or negative). The physical manifestation of this topological quantization is that the magnetic skyrmions are topologically protected from decaying to a uniform magnetization that has a winding number equal to zero. As we will show in the sections dedicated to results and discussion, the spins from the core of the skyrmion are protected by an energy barrier driven by the DMI from the reversed spins from the edge.

Ab-Initio Calculations of Magnetic Properties
The ab-initio calculations for extracting and analyzing the magnetic properties have been performed using the Full Potential Linear Augmented Plane Wave FP-LAPW Wien2k code [23]. The ab-initio results have been used, in a first step, as starting points for further micromagnetic calculations. For modeling multilayer stacks compatible with PMA and DMI sizeable effects we used supercells with typical X/Fe(t Fe )/MgO(001) configuration, where X = nonmagnetic metal (Au, V, Ag) or insulator (MgO). The thickness of the ferromagnetic Fe layer (t Fe ) was chosen in a few monolayers (ML) range (e.g., 5ML), as required for promoting the perpendicular magnetization ground states driven by the interplay of both top and bottom interfacial PMA. The calculation of the magnetic anisotropy energy (MAE) was executed by means of a fully relativistic spin-orbit scheme using the total energy and force theorem approaches [24], both providing similar results. Using the total energy approach, including the spin-orbit, the perpendicular magnetic anisotropy (PMA) has been estimated as the total energy difference between the easy (magnetization perpendicular to the film plane) and the hard directions (in-plane magnetization). Considering the extreme sensitivity of the magnetic anisotropy energy to the k-space meshing, first, a convergence study of the total energy with respect to the total number of k-points has been thoroughly performed.
The spin-orbit term of the non-relativistic Dirac (Pauli) Hamiltonian: ∂V ∂z is the Rashba constant which is a measure of the spin-orbit interaction and → e z the unit vector of the Oz (which corresponds to the electron confinement direction). According to its definition, one can conclude that α R is proportional to the electric field E = − ∂V ∂z pointing along the Oz direction. In a multilayer stack, between two adjacent materials, such intrinsic electric fields naturally exist and lead to Rashba spinorbit interaction effects. This bidimensional case is a good description of ultrathin films, their thickness in the few monolayer range being smaller than the quantum coherence length of electrons. In the situation when a metal-insulator (or metal semiconductor) are associated in a multilayer heterostructure, forming an interface, a depletion zone is defined with a significant interfacial electric field associated. By diagonalizing the Rashba Hamiltonian, the eigenvalues will be: E ± (k I I ) = 2 k 2 I I 2m 0 ± α R |k I I |, which correspond to parabolic bands with an offset of the parabola minimum in positive (or negative) k values (Rashba splitting). The minimum of the parabola can be found when derivative condition is Considering this expression of E 0 , we concluded that α R could be calculated using the values of E and k corresponding to the minimum of the parabolic dispersion band: α R = 2E 0 k 0 . We used this expression to deduce the Rashba parameter from the Rashba parabolic offset of the bands calculated for different supercell architectures describing magnetic heterostructures. Knowing the Rashba parameter, one can calculate the Rashba contribution to the DMI as [25]: DMI = 2k R A, with A = the exchange stiffness of the ferromagnetic material (e.g., 21 pJ/m for Fe), k R = 2α R m e 2 , with m e the electron effective mass of the Rashba shifted parabolic bands.
Within these formalisms, the electric field (E-field) can be applied using a zig-zag additional potential in the Hamiltonian (Figure 3), as implemented by Stahn et al. [26]. Therefore, for different values of the electric field, the (perpendicular) magnetic anisotropy (PMA), and the Dzyaloshinskii-Moriya interaction (DMI) can be calculated to provide theoretical insight on their origin and electric field dependence.
The effective field can be calculated as the functional derivative of the magnetic free energy of the system E[m,t] that includes various contributions: Zeemann, demagnetizing, anisotropy, exchange: symmetric (J0) and asymmetric (DMI), STT effects, etc. : = − ð ð . In our calculations, we have employed an NVIDIA Geforce RTX 3070 graphic card with 5888 CUDA cores for parallel computing. The modeled systems are nanometric-size disks constituted from tri-layer stacks, like those modeled by ab-initio calculations, in which a ferromagnetic film is sandwiched between a top (T) and a bottom (B) nonmagnetic layer. Therefore, the PMA and DMI are controlled by additive interface effects (top panel of Figure 3). In the Mumax 3 code, the perpendicular magnetic anisotropy was introduced as a 1st order uniaxial anisotropy Ku along the Oz (001) axis, and the DMI has been chosen to have interfacial origin (Dind), as expected for ultrathin films stacks. In the simulations, we included as variables: other specific material/multilayer stack parameters such as the α-Gilbert damping of FM, exchange stiffness Aex, the saturation magnetization of the FM, the shape/size of the patterned nanostructure (e.g., disk, track, …), the density of a spin-polarized current injected in a nano-disk for STT writing of skyrmions and the time Therefore, for different values of the electric field, the (perpendicular) magnetic anisotropy (PMA), and the Dzyaloshinskii-Moriya interaction (DMI) can be calculated to provide theoretical insight on their origin and electric field dependence.

Micromagnetic Modeling Tools
The micromagnetic simulations have been performed using the Mumax 3 GPU accelerated code (Ghent University, Ghent, Belgium) [28], based on solving the Landau-Lifshitz-Gilbert equation: where B e f f is the effective field, α is the Gilbert damping parameter, γ is the gyromagnetic ratio (1.75 × 10 7 s Oe −1 ) and µ 0 the vacuum permeability (4π × 10 −7 F/m). The effective field can be calculated as the functional derivative of the magnetic free energy of the system E[m,t] that includes various contributions: Zeemann, demagnetizing, anisotropy, exchange: symmetric (J 0 ) and asymmetric (DMI), STT effects, etc. : In our calculations, we have employed an NVIDIA Geforce RTX 3070 graphic card with 5888 CUDA cores for parallel computing. The modeled systems are nanometric-size disks constituted from tri-layer stacks, like those modeled by ab-initio calculations, in which a ferromagnetic film is sandwiched between a top (T) and a bottom (B) nonmagnetic layer. Therefore, the PMA and DMI are controlled by additive interface effects (top panel of Figure 3). In the Mumax 3 code, the perpendicular magnetic anisotropy was introduced as a 1st order uniaxial anisotropy K u along the Oz (001) axis, and the DMI has been chosen to have interfacial origin (D ind ), as expected for ultrathin films stacks. In the simulations, we included as variables: other specific material/multilayer stack parameters such as the α-Gilbert damping of FM, exchange stiffness A ex , the saturation magnetization of the Nanomaterials 2022, 12, 4411 8 of 33 FM, the shape/size of the patterned nanostructure (e.g., disk, track, . . . ), the density of a spin-polarized current injected in a nano-disk for STT writing of skyrmions and the time length of current or gating electric field pulses. For dynamical simulations, the LLG equation has been integrated using the adaptive step RK45, the Dormand-Prince solver. For the calculations of the phase diagrams, the energy minimum of the states has been obtained using the mumax3 relax () function that disables the precession term of the LLG equation allowing the system to relax towards the minimum of energy. For including the temperature effects, a finite temperature has been set to the ground state after the relaxation. Then, the evolution/stability of the system has been studied within a finite time integration window. The classification of the final micromagnetic states has been performed by topological charge calculations: dy tested against either artificial intelligence (AI) algorithms (image recognition).

Ab-initio Calculations of Multilayer Stacks
These calculations provide direct insight into the magnetic properties responsible for the skyrmion stabilization and their control by external stimuli, e.g., electric or magnetic fields. These parameters are the perpendicular magnetic anisotropy PMA, the asymmetric exchange DMI and the saturation magnetization, calculated as the total magnetic moment divided by the volume of the unit supercell.
First, we focused on the charge depletion effects at the interfaces between the ferromagnetic film and the bottom X and top MgO layers. Figure 4 illustrates a typical example of charge density representation, calculated for an Au/Fe/MgO supercell. A similar analysis was performed for other bottom underlayers: e.g., X = V, Ag, MgO. Corresponding to this charge depletion, an intrinsic electric field would appear at each interface. Its orientation depends on the relative magnitude of the work functions in Fe and X and Fe and MgO, respectively. In the case of the Au/Fe/MgO stack, from the valence charge density distribution one can distinguish the intrinsic charge depletion zones at the bottom Au/Fe and top Fe/MgO interfaces, and, from the color map corresponding to the charge density variation ∆n, the orientation and the magnitude of the corresponding intrinsic electric fields E B and E T , that would be responsible for the PMA and DMI. Moreover, one can qualitatively observe that a positive (negative) voltage bias has an effect, the decrease (enhancement) of the intrinsic field at the top Fe/MgO interface ( Figure 4c) and, as expected, no effect for the metallic Au/Fe bottom interface ( Figure 4d). It was demonstrated that this intrinsic interface electric field would be responsible for both PMA and DMI [25,29]. Moreover, as theoretically predicted by Barnes et Maekawa [30], the additive effect of the intrinsic electric fields of both top and bottom interfaces would lead to a net Rashba contribution to the PMA and DMI whose magnitude and signs can be modulated by the chemical nature of the underlayer and overlayer. Therefore, within this interface Rashba framework, one can easily understand the effect of an external gating electric field which, depending on its orientation will either enhance or decrease the value of the intrinsic interface fields and implicitly, the magnitude of the PMA and the DMI.
Indeed, our ab-initio calculations fully subscribe to this framework and confirm the quantum analytical predictions of Barnes et al. [30], illustrating the opposite sign variation of the PMA with the electric field, corresponding to V/Fe(5ML)/MgO and Au/Fe(5ML)/MgO systems ( Figure 5). The PMA, described by the effective anisotropy constant K u (J/m 3 ), is calculated from the surface anisotropy energy K s (J/m 2 ), using the equation K u (J/m 3 ) = K s (J/m 2 )/ t Fe , with t Fe being the thickness of Fe (0.715 nm for 5ML). The value of K s is obtained by dividing the total energy difference between the easy (magnetization perpendicular to the film plane) and the hard directions (in-plane magnetization) to the area of the supercell base ( Figure 4a). Indeed, our ab-initio calculations fully subscribe to this framework and confirm the quantum analytical predictions of Barnes et al. [30], illustrating the opposite sign variation of the PMA with the electric field, corresponding to V/Fe(5ML)/MgO and Au/Fe(5ML)/MgO systems ( Figure 5). The PMA, described by the effective anisotropy constant Ku (J/m 3 ), is calculated from the surface anisotropy energy Ks (J/m 2 ), using the equation Ku (J/m 3 ) = Ks (J/m 2 )/ tFe, with tFe being the thickness of Fe (0.715 nm for 5ML). The value of Ks is obtained by dividing the total energy difference between the easy (magnetization perpendicular to the film plane) and the hard directions (in-plane magnetization) to the area of the supercell base ( Figure 4a). This opposite behavior can be explained by considering the effect on an external electric field on the additive Rashba field of each interface in the X/Fe/MgO stack. In Figure 6 we illustrate the potential profiles seen by an electron in the X/Fe/MgO stack, considering the corresponding work functions of the three elements: X (V and Au), Fe and MgO. These profiles agree with the orientation of the intrinsic electric fields at each interface ⃗ , i = 1,2 obtained by analyzing the interfacial depletion of the electronic charge density, similar to This opposite behavior can be explained by considering the effect on an external electric field on the additive Rashba field of each interface in the X/Fe/MgO stack. In Figure 6 we illustrate the potential profiles seen by an electron in the X/Fe/MgO stack, considering the corresponding work functions of the three elements: X (V and Au), Fe and MgO. These profiles agree with the orientation of the intrinsic electric fields at each interface → E i , i = 1,2 obtained by analyzing the interfacial depletion of the electronic charge density, similar to the situation shown in Figure 4. This opposite behavior can be explained by considering the effect on an external electric field on the additive Rashba field of each interface in the X/Fe/MgO stack. In Figure 6 we illustrate the potential profiles seen by an electron in the X/Fe/MgO stack, considering the corresponding work functions of the three elements: X (V and Au), Fe and MgO. These profiles agree with the orientation of the intrinsic electric fields at each interface ⃗ , i = 1,2 obtained by analyzing the interfacial depletion of the electronic charge density, similar to the situation shown in Figure 4. The corresponding interfacial Rashba magnetic fields are ⃗ = − ( ⃗ × ⃗ ), the net (additive) Rashba field felt by the Fe layer, ⃗ = ⃗ + ⃗ would have a magnitude whose sign would depend on the type of the X underlayer (as confirmed by band structure analysis depicted in Figure 7). Therefore, the effect of an external electric field would be different. Depending on the nature of X, a net increase/decrease in the ⃗ (and, consequently, of the Rashba contribution to the PMA and DMI) would be observed. In the case of Au/Fe/MgO, from the valence charge density analysis ( Figure 4c) we saw that a positive bias voltage (corresponding to a negative external electric field) leads to a reduction in the top interface electric field and, consequently, to a reduction in MAE, as seen in Figure 5b).
Furthermore, the DMI can be calculated from the Rashba parameter , estimated from the shift in k of the band structure corresponding to a spin-orbit ab-initio calculation with the magnetization M fixed along (100) and (−100) directions ( Figure 7). The corresponding interfacial Rashba magnetic fields are net (additive) Rashba field felt by the Fe layer, would have a magnitude whose sign would depend on the type of the X underlayer (as confirmed by band structure analysis depicted in Figure 7). Therefore, the effect of an external electric field would be different. Depending on the nature of X, a net increase/decrease in the → B net R (and, consequently, of the Rashba contribution to the PMA and DMI) would be observed. In the case of Au/Fe/MgO, from the valence charge density analysis ( Figure 4c) we saw that a positive bias voltage (corresponding to a negative external electric field) leads to a reduction in the top interface electric field and, consequently, to a reduction in MAE, as seen in Figure 5b). Considering this framework, our scope was to compute, explain and correlate the sign, the magnitude and the electric field dependence of the anisotropy, the Rashba coefficient , and the Dzyaloshinskii-Moriya (DMI) interaction parameter in different types of X/Fe/MgO(001) systems in which the variable layer was X (where X = V, Au, Ag, Pd, MgO, etc.). A typical example of such kind of result is depicted in Figure 8, corresponding Furthermore, the DMI can be calculated from the Rashba parameter α R , estimated from the shift in k of the band structure corresponding to a spin-orbit ab-initio calculation with the magnetization M fixed along (100) and (−100) directions ( Figure 7).
Considering this framework, our scope was to compute, explain and correlate the sign, the magnitude and the electric field dependence of the anisotropy, the Rashba coefficient α R , and the Dzyaloshinskii-Moriya (DMI) interaction parameter in different types of X/Fe/MgO(001) systems in which the variable layer was X (where X = V, Au, Ag, Pd, MgO, etc.). A typical example of such kind of result is depicted in Figure 8, corresponding to the ab-initio analysis of the Au/Fe(5ML)/MgO system. in the two situations (adapted from [28]).
Considering this framework, our scope was to compute, explain and correlate the sign, the magnitude and the electric field dependence of the anisotropy, the Rashba coefficient , and the Dzyaloshinskii-Moriya (DMI) interaction parameter in different types of X/Fe/MgO(001) systems in which the variable layer was X (where X = V, Au, Ag, Pd, MgO, etc.). A typical example of such kind of result is depicted in Figure 8, corresponding to the ab-initio analysis of the Au/Fe(5ML)/MgO system. We also studied the influence of a monolayer of different metals, e.g., Pt, Au or Pd inserted at the top interface between Fe and MgO. We observed that in the case of a Pt, the PMA, the DMI and their variation with respect to an external electric field are always enhanced, independently of the nature of X. The most significant effect has been found for the Au/Fe(5ML)/Pt(1ML)/MgO system. Here, the Pt ad-monolayer amplifies the PMA by an order of magnitude, from 1.91 to 15.77 MJ/m 3 , the DMI from 1.67 to 2.27 mJ/m 2 , and the variation of PMA with E-field roughly by a factor of 2. These results indicate promising strategies for tailoring the PMA, DMI and their response to an electric field, very important issues for micromagnetic properties and control, as we will illustrate later.
Therefore, the ab-initio calculations will be very important for a further micromagnetic analysis of skyrmionic nanomaterials. The opposite sign of the Rashba coefficient in different systems would determine the opposite signs of the Dzyaloshinskii-Moriya Figure 8. Electric field dependence of the surface anisotropy K u (PMA), α R and DMI, here for the Au/Fe(5ML)/MgO system-(adapted from [28]). The zero field values: K u (E-field = 0) = 1.91 MJ/m 3 and DMI (E-field = 0) = 1.67 mJ/m 2 will be further used in micromagnetic simulations (see Section 3.2.2).
We also studied the influence of a monolayer of different metals, e.g., Pt, Au or Pd inserted at the top interface between Fe and MgO. We observed that in the case of a Pt, the PMA, the DMI and their variation with respect to an external electric field are always enhanced, independently of the nature of X. The most significant effect has been found for the Au/Fe(5ML)/Pt(1ML)/MgO system. Here, the Pt ad-monolayer amplifies the PMA by an order of magnitude, from 1.91 to 15.77 MJ/m 3 , the DMI from 1.67 to 2.27 mJ/m 2 , and the variation of PMA with E-field roughly by a factor of 2. These results indicate promising strategies for tailoring the PMA, DMI and their response to an electric field, very important issues for micromagnetic properties and control, as we will illustrate later.
Therefore, the ab-initio calculations will be very important for a further micromagnetic analysis of skyrmionic nanomaterials. The opposite sign of the Rashba coefficient α R in different systems would determine the opposite signs of the Dzyaloshinskii-Moriya interaction (DMI). This is particularly important because the sign of the DMI determines the chirality of the chiral structures (domain-walls, skyrmions) and, in correlation with the orientation of the core polarization, their velocity direction driven by spin currents (STT). Moreover, the sign and the magnitude of the PMA and DMI variation with an electric field determines the possible trajectory for the displacement within the magnetic parameters phase diagrams (discussed in the next section) under the applied electric field, in view of skyrmions stabilization and manipulation for various applications.

Micromagnetic Design of Skyrmionic Materials
This section addresses the question about the critical material magnetic parameters for injecting, hosting as ground states, and manipulation of skyrmions in patterned nanostructures. Using micromagnetic simulations, we calculated different types of phase diagrams of final magnetic states in nanometric-size patterned disks. First, we used as input the magnetic parameters obtained from ab-initio calculations for the Au/Fe/MgO system, in which the PMA, DMI and their variation with gating electric fields are controlled by interfacial cumulative effects (Figure 6b). Within these circumstances, we calculated the phase diagrams for skyrmion injection by spin-transfer-torque (STT) effects related to spinpolarized currents in patterned nanodisks of various diameters. Second, the skyrmionic injection by STT has been addressed via extended phase diagrams, in which the current pulse parameters have been fixed and the material parameters, PMA and DMI, have been varied in a wider range, to cover other potential skyrmionic materials. Then, to get a deeper understanding of mechanisms and magnetic free-energy contributions responsible for the stabilization and the stability of skyrmions, our study has been further extended to phase diagrams of ground states, in nanodisks with different sizes and shapes. We succeeded to identify the occurrence of various micromagnetic chiral states and analyzed the skyrmions stability/evolution as a function of magnetic and geometric parameters: PMA, DMI, M s and disk size. We briefly extrapolated the phase diagrams to the case of antiferromagnetically coupled skyrmions that could be stabilized in synthetic antiferromagnetic materials with interesting perspectives for racetrack memories and coupled qubits applications. Furthermore, we discussed some skyrmion manipulation issues: in nano-oscillators and racetracks, we illustrate the possibility to generate multi-skyrmionic ground states interesting for memristor functions in neuromorphic devices. Finally, we address some issues related to the skyrmionic state manipulation (core polarization and chirality by external gating electric fields). The choice of the diameter size of the disk in most simulations (90 nm) has been optimized as a compromise between the following issues: (i) our initial fixed target, to have nanometric size skyrmions in materials with experimentally achievable DMI (having in view our observation concerning the inverse proportionality between the size of a skyrmion and the DMI required to stabilize it), (ii) the efficiency and versatility of the skyrmionic injection by STT in larger disks increases with the disk size.

Injection of Skyrmions by STT in Patterned Nanodisks
In this paragraph, we address the possibility of stabilizing skyrmions in patterned disks by writing them from a saturated ground state, using the STT effect of a spin-polarized current injected perpendicularly across the disk. The geometry used for this micromagnetic simulation is represented in Figure 9. It is compatible with a magnetic tunnel junction (MTJ) element, either the individual or as the writing/reading part in a racetrack memory device. Within the hard-soft architecture of the MTJ, the pillar structure has the following parts: (1) a bottom ferromagnetic disk, constituted from the soft magnetic material in which the skyrmions will be written (free magnetic layer), (2) an insulating tunnel barrier, and (3) a top fixed (hard) magnetic layer acting as a spin polarizer for the current injected across the insulator. For the Mumax3 micromagnetic simulations we used the following geometry: 1 nm thick circular disks of 90 nm diameter for the free (bottom) layer and 20 nm diameter for the hard (top) polarizer, respectively. The fixed and spacer layer parameters in Mumax3 code, required for including the Slonczewski and Zhang-Li spin transfer torques effects [27] were set up as follows: the Slonczewski Λ parameter λ = 1, the electric current polarization Pol = 0.4-1, the Slonczewski secondary STT term ε' = 0, the non-adiabaticity parameter of the Zhang-Li model ξ = 2, and the magnetization of the fixed layer was locked along the (001) perpendicular direction. The material parameters for the free layer in which the skyrmions will be written were chosen to be exactly those issued from the ab-initio calculations for the Au/Fe(5ML)/MgO system: K u (PMA) = 1.91 × 10 +6 J/m 3 and D ind (DMI) = 1.67 × 10 −3 J/m 2 ; M s = 1714 × 10 +3 A/m and α = 0.01 corresponding to the bulk Fe. In our simulations, we also considered as a variable parameter the voltage modulation of the anisotropy ( ξ VCMA ) and the DMI for the bottom free layer. During the current pulse, the MTJ gets polarized and the electric field across the insulator was chosen to reduce the PMA and the DMI: K u = Ku ini (1-ξ VCMA ) and D ind = D ind ini (1-ξ VCMA ); the initial Ku ini and DMI ini values being restored after the pulse. The size of the unit cell of the micromagnetic grid was 1 × 1 × 1 nm. 10 J/m and Dind (DMI) = 1.67 × 10 J/m ; Ms = 1714 × 10 A/m and α = 0.01 corresponding to the bulk Fe. In our simulations, we also considered as a variable parameter the voltage modulation of the anisotropy ( ξVCMA) and the DMI for the bottom free layer. During the current pulse, the MTJ gets polarized and the electric field across the insulator was chosen to reduce the PMA and the DMI: Ku = Ku ini (1-ξVCMA) and Dind = Dind ini (1-ξVCMA); the initial Ku ini and DMI ini values being restored after the pulse. The size of the unit cell of the micromagnetic grid was 1 × 1 × 1 nm. Within this framework, we have calculated various Tpulse-Jc-Q and Tpulse-Jc-mz phase diagrams, starting either from an initial antiparallel (AP) or parallel (P) state of the fixed polarizer and the free layer; Tpulse is the length of the current pulse, Jc the density of the polarized electric current, Q the topological charge and mz the normalized z component of the magnetization. The micromagnetic grid used is a unit cell of 1 × 1 × 1 nm and the diagrams are built from 10.000 relaxed dynamics. In Figure 10, we illustrate the skyrmion writing phase diagrams calculated considering a current polarization p = 0.8 and a voltage modulation of anisotropy and DMI, ξ = 0.1. Like in common experiments of magnetization switch by STT, as a function of the initial state, P or AP, the direction of the injected Jc for writing the skyrmions by STT must be reversed (see Figure 10a,d). The Tpulse-Jc-Q diagrams (Figure 10b,d), corresponding to injection from AP(P) initial states, illustrate the Jc and Tpulse parameter range allowing to write a skyrmion: either with a positive core polarization (Figure 10a), when STT writing from an initial AP saturated state or a negative core polarization (Figure 10d), when injecting from an initial P state. Analyzing the sign and the value of the topological charge of the skyrmionic states range of the Tpulse-Jc-Q diagrams: (Q = 1, red color in (b) or Q = -1, blue color in (e)), one can see that, depending on the initial AP or P state, the core polarization, and the chirality of the written skyrmion will be opposite. Moreover, the Tpulse-Jc-mz phase diagrams (Figure 10c,d) bring complementary information, indistinguishable from the Tpulse-Jc-Q diagrams. They allow us to discriminate the parameter range Jc -Tpulse leading to the full reversal of the free-layer disk by STT, the STT writing experiments in STT-RAMs are the same. Within this framework, we have calculated various T pulse -J c -Q and T pulse -J c -m z phase diagrams, starting either from an initial antiparallel (AP) or parallel (P) state of the fixed polarizer and the free layer; T pulse is the length of the current pulse, J c the density of the polarized electric current, Q the topological charge and m z the normalized z component of the magnetization. The micromagnetic grid used is a unit cell of 1 × 1 × 1 nm and the diagrams are built from 10.000 relaxed dynamics. In Figure 10, we illustrate the skyrmion writing phase diagrams calculated considering a current polarization p = 0.8 and a voltage modulation of anisotropy and DMI, ξ = 0.1. Like in common experiments of magnetization switch by STT, as a function of the initial state, P or AP, the direction of the injected J c for writing the skyrmions by STT must be reversed (see Figure 10a,d). The T pulse -J c -Q diagrams (Figure 10b,d), corresponding to injection from AP(P) initial states, illustrate the J c and T pulse parameter range allowing to write a skyrmion: either with a positive core polarization (Figure 10a), when STT writing from an initial AP saturated state or a negative core polarization (Figure 10d), when injecting from an initial P state. Analyzing the sign and the value of the topological charge of the skyrmionic states range of the T pulse -J c -Q diagrams: (Q = 1, red color in (b) or Q = -1, blue color in (e)), one can see that, depending on the initial AP or P state, the core polarization, and the chirality of the written skyrmion will be opposite. Moreover, the T pulse -J c -m z phase diagrams (Figure 10c,d) bring complementary information, indistinguishable from the T pulse -J c -Q diagrams. They allow us to discriminate the parameter range J c -T pulse leading to the full reversal of the free-layer disk by STT, the STT writing experiments in STT-RAMs are the same.
Similar phase diagrams, calculated for disks with various sizes, illustrate that the efficiency of the STT writing increases with increasing the disk diameter (see Figure 11). The J c and T pulse window for the (S) writing by STT is significantly larger for the 90 nm disks as compared to 60 nm ones.
Based on the information brought by the T pulse -J c -Q and T pulse -J c -m z phase diagrams, we have chosen a specific point whose parameters (J c , T pulse ) are suitable for writing a skyrmionic state by STT (white point designed by (WP) in Figure 10b,e). Corresponding to these (J c , T pulse ) values that were kept fixed, we have calculated new phase diagrams in which the PMA and the DMI have been varied to potentially cover a broader range corresponding to other skyrmionic materials. The diagrams illustrated in Figure 12, correspond to a 90 nm diameter circular nanodisk, patterned from a 1 nm thick ferromagnetic film, and are built from 10.000 relaxed dynamics in which we chose as variables the anisotropy and the asymmetric exchange (DMI) and fixed the other magnetic parameters: the saturation magnetization M s = 1714 kA/m; the exchange stiffness A ex = 2.1 × 10 +11 J/m; the Gilbert damping constant α = 0.01, which are typical for a Fe thin film that can enter in a X/Fe(5ML)/MgO multilayer stack. The color scale represents the value of the topological used to classify the micromagnetic configuration of the final relaxed chiral ground states. In our simulations, we start from an initial saturated P or AP (with respect to the hard reference layer) state, and after a current pulse (J c , T pulse ) we cut the pulse and determine the final relaxed micromagnetic configuration of the disk. Similar phase diagrams, calculated for disks with various sizes, illustrate that the efficiency of the STT writing increases with increasing the disk diameter (see Figure 11). The Jc and Tpulse window for the (S) writing by STT is significantly larger for the 90 nm disks as compared to 60 nm ones. Based on the information brought by the Tpulse-Jc-Q and Tpulse-Jc-mz phase diagrams, we have chosen a specific point whose parameters (Jc, Tpulse) are suitable for writing a skyrmionic state by STT (white point designed by (WP) in Figure 10b,e). Corresponding to these (Jc, Tpulse) values that were kept fixed, we have calculated new phase diagrams in which the PMA and the DMI have been varied to potentially cover a broader range corresponding to other skyrmionic materials. The diagrams illustrated in Figure 12, correspond to a 90 nm diameter circular nanodisk, patterned from a 1 nm thick ferromagnetic film, and are built from 10.000 relaxed dynamics in which we chose as variables the anisotropy and the asymmetric exchange (DMI) and fixed the other magnetic parameters: the satura- Similar phase diagrams, calculated for disks with various sizes, illustrate that the efficiency of the STT writing increases with increasing the disk diameter (see Figure 11). The Jc and Tpulse window for the (S) writing by STT is significantly larger for the 90 nm disks as compared to 60 nm ones. Based on the information brought by the Tpulse-Jc-Q and Tpulse-Jc-mz phase diagrams, we have chosen a specific point whose parameters (Jc, Tpulse) are suitable for writing a skyrmionic state by STT (white point designed by (WP) in Figure 10b,e). Corresponding to these (Jc, Tpulse) values that were kept fixed, we have calculated new phase diagrams in which the PMA and the DMI have been varied to potentially cover a broader range corresponding to other skyrmionic materials. The diagrams illustrated in Figure 12, correspond to a 90 nm diameter circular nanodisk, patterned from a 1 nm thick ferromagnetic film, and are built from 10.000 relaxed dynamics in which we chose as variables the anisotropy and the asymmetric exchange (DMI) and fixed the other magnetic parameters: the saturation magnetization Ms = 1714 kA/m; the exchange stiffness Aex = 2.1 × 10 +11 J/m; the Gilbert These diagrams, illustrated in Figure 12, show that the stabilization of a skyrmionic state occurs in an optimum window of DMI; the larger the DMI, the wider the PMA window will be in which skyrmionic states can be injected by a current. Above a certain DMI value, the skyrmionic states become unstable in a PMA range increasing with increasing DMI, and other complex chiral structures will be stabilized. When writing the skyrmion by STT using a current pulse, above a certain DMI value the skyrmionic final states become also unfavorable. After the current pulse, the system will relax towards a vortex, a perpendicularly saturated or a complex chiral final state. This will happen within a PMA window whose width increases with the increase in the DMI magnitude. damping constant α = 0.01, which are typical for a Fe thin film that can enter in a X/Fe(5ML)/MgO multilayer stack. The color scale represents the value of the topological charge: = ⃗ • ⃗ × ⃗ used to classify the micromagnetic configuration of the final relaxed chiral ground states. In our simulations, we start from an initial saturated P or AP (with respect to the hard reference layer) state, and after a current pulse (Jc, Tpulse) we cut the pulse and determine the final relaxed micromagnetic configuration of the disk. These diagrams, illustrated in Figure 12, show that the stabilization of a skyrmionic state occurs in an optimum window of DMI; the larger the DMI, the wider the PMA window will be in which skyrmionic states can be injected by a current. Above a certain DMI value, the skyrmionic states become unstable in a PMA range increasing with increasing DMI, and other complex chiral structures will be stabilized. When writing the skyrmion by STT using a current pulse, above a certain DMI value the skyrmionic final states become also unfavorable. After the current pulse, the system will relax towards a vortex, a perpendicularly saturated or a complex chiral final state. This will happen within a PMA window whose width increases with the increase in the DMI magnitude.
To get a deeper understanding of the mechanisms and magnetic free-energy contributions responsible on the stabilization and the stability of skyrmionic states, we extended our study to phase diagrams of ground states in nanodisks with different sizes and shapes.

Skyrmionic Ground States in Patterned Nanodisks
The ground state of the magnetic nanostructure has been obtained by relaxing the LLG magnetization dynamics, starting from an initial saturated state, chosen with the magnetization aligned along the -Oz = (00-1) direction. Figure 13 illustrates such a diagram, corresponding to a 90 nm diameter circular nanodisk, patterned from a 1 nm thick ferromagnetic film. The parameters used in these simulations are the same as those used for generating the diagram from Figure 12. Starting from this diagram, we analyzed the evolution of the ground state micromagnetic configurations when fixing a parameter and varying the other (e.g., PMA = 1.5 × 10 +6 J/m 3 , DMI variable-line (1)) and DMI = 3 × 10 −3 J/m 2 , variable PMA-line (2). The resulting magnetic features are displayed in Figures 13  and 14 and the following conclusions have been drawn: Figure 12. DMI-PMA-Q phase diagrams corresponding to simulation experiments started either from an initial P state (a) or an initial AP state (b). The parameters of the current pulse used in the simulation were: T pulse = 0.5 ns and J c = −3 × 10 12 A/m 2 or J c = + 3 × 10 12 A/m 2 , depending on the initial P or AP state (WP in Figure 22a,d). The sign of the written skyrmion topological charge is opposite in the two situations.
To get a deeper understanding of the mechanisms and magnetic free-energy contributions responsible on the stabilization and the stability of skyrmionic states, we extended our study to phase diagrams of ground states in nanodisks with different sizes and shapes.

Skyrmionic Ground States in Patterned Nanodisks
The ground state of the magnetic nanostructure has been obtained by relaxing the LLG magnetization dynamics, starting from an initial saturated state, chosen with the magnetization aligned along the -Oz = (00-1) direction. Figure 13 illustrates such a diagram, corresponding to a 90 nm diameter circular nanodisk, patterned from a 1 nm thick ferromagnetic film. The parameters used in these simulations are the same as those used for generating the diagram from Figure 12. Starting from this diagram, we analyzed the evolution of the ground state micromagnetic configurations when fixing a parameter and varying the other (e.g., PMA = 1.5 × 10 +6 J/m 3 , DMI variable-line (1)) and DMI = 3 × 10 −3 J/m 2 , variable PMA-line (2). The resulting magnetic features are displayed in Figures 13 and 14 and the following conclusions have been drawn: • At a fixed DMI, at small PMA values, the ground state is vortex-like, with Q ≈ −0.5; when the PMA increases, within a certain range whose extension increases with the value of the DMI, the ground state becomes skyrmionic Q ≈ −1; above a threshold maximum value of the PMA, the ground state will be a saturated monodomain.

•
Within the skyrmionic window, the skyrmion's diameter progressively decreases with the PMA increase. The narrow transition from a skyrminonic to a saturated state implies some complex chiral states with |Q| > 1.

•
The skyrmionic states are stabilized by larger values of DMI but, above a critical value, at fixed PMA, large DMI will favour |Q| > 1 configurations within a PMA window increasing with the increase in the DMI.

•
In conclusion, an optimum window of PMA and DMI values is required to stabilize a skyrmionic ground state.
The magnetization vector field projection on the Bloch unit sphere of the Néel chiral structures ("hedgehog") illustrates the correlation between the value of the topological charge Q and the "population" of the Bloch sphere by spins. For a vortex state, half of the Bloch sphere is populated, a skyrmionic magnetization vector field populates the entire sphere and a saturated state populates either the (0;0;−1) or the (0;0;1) poles of the unit sphere m 2 x + m 2 y + m 2 z = 1. • At a fixed DMI, at small PMA values, the ground state is vortex-like, with Q ≈ −0.5; when the PMA increases, within a certain range whose extension increases with the value of the DMI, the ground state becomes skyrmionic Q ≈ −1; above a threshold maximum value of the PMA, the ground state will be a saturated monodomain.

•
Within the skyrmionic window, the skyrmion's diameter progressively decreases with the PMA increase. The narrow transition from a skyrminonic to a saturated state implies some complex chiral states with | | > 1.

•
The skyrmionic states are stabilized by larger values of DMI but, above a critical value, at fixed PMA, large DMI will favour | | > 1 configurations within a PMA window increasing with the increase in the DMI.

•
In conclusion, an optimum window of PMA and DMI values is required to stabilize a skyrmionic ground state.
The magnetization vector field projection on the Bloch unit sphere of the Néel chiral structures ("hedgehog") illustrates the correlation between the value of the topological charge Q and the "population" of the Bloch sphere by spins. For a vortex state, half of the Bloch sphere is populated, a skyrmionic magnetization vector field populates the entire sphere and a saturated state populates either the (0;0;−1) or the (0;0;1) poles of the unit sphere + + = 1.
The total energy of the final states is different. The analysis of a line profile in the total magnetic energy density of the disk demonstrates that the vortex state is less energetically favorable for the inner spins and that the skyrmionic state is a topologically stabilized configuration (Figure 15c). Here, for the spins in the inner skyrmionic core, an energy barrier ∆ prevents them to rotate towards the opposite orientation of the edge spins. The height of this barrier is a complex energetical competition issue between the favorable DMI and PMA contributions and unfavorable magnetostatic and direct exchange. For room temperature skyrmion stability, the barrier energy has to be larger than the thermal energy , being the Boltzmann constant and T the value of the absolute temperature (more details can be The total energy of the final states is different. The analysis of a line profile in the total magnetic energy density of the disk demonstrates that the vortex state is less energetically favorable for the inner spins and that the skyrmionic state is a topologically stabilized configuration (Figure 15c). The variation of the energy density barrier ∆ separating the skyrmion core spins and the oppositely oriented edge ones has a non-monotonous evolution with DMI and PMA, as illustrated in Figure 16. In both situations, an optimum window exists for stable skyrmionic states (Figure 16b,c): complex chiral structures with Q > 0 or saturated state (Q = 0) appearing after some up-threshold DMI and PMA values, respectively. A larger PMA implicates a confinement of the skyrmion core in a smaller region. This would require a larger DMI for twisting the spins against the direct exchange interaction. The threshold from (S) to (P) is roughly showing a parabolic increase with the DMI (Figure  15a) dotted white line.
A deeper insight into the skyrmion stability, with respect to the increase in the anisotropy, can be obtained by analyzing the corresponding variation of different contributions to the magnetic free energy (Figure 17). The analysis has been performed along a path, indicated by a solid black vertical line in Figure 15, corresponding to a fixed value of DMI (Dind = 4 mJ/m 2 ) and variable PMA (Ku = 0− 2.5 × 10 +6 J/m 3 ). We see that the skyrmion stability is a complex issue, being determined by the balance between the PMA, the DMI, the direct exchange and the demagnetizing (magnetostatic) energies. When the PMA increases during the evolution from a vortex state to a skyrmionic one, the total energy decreases. However, this evolution costs more and more in demagnetizing energy which is continuously increasing. At a certain critical value of PMA, an abrupt transition into a saturated perpendicular state will stop the increase in the magnetostatic energy and will accommodate both the anisotropy and the direct exchange energies, in detriment to the DMI energy (the abrupt transition from (S) to (P) is also illustrated in Figure 17). Here, for the spins in the inner skyrmionic core, an energy barrier ∆E dens prevents them to rotate towards the opposite orientation of the edge spins. The height of this barrier is a complex energetical competition issue between the favorable DMI and PMA contributions and unfavorable magnetostatic and direct exchange. For room temperature skyrmion stability, the barrier energy has to be larger than the thermal energy k B T, k B being the Boltzmann constant and T the value of the absolute temperature (more details can be found in Appendix A).
Within the skyrmionic range of the phase diagram ( Figure 15) the increase in either PMA or DMI will lead to the shrinkage of the skyrmionic core.
When the DMI increases at fixed PMI (path (1)), the skyrmonic size decreases, but above a critical DMI threshold complex chiral structures with Q > 1 become stable. In the case of the PMA increase at fixed DMI (path (2)), above a critical threshold, the core shrinks to collapse into a saturated state with an opposite magnetization orientation with respect to the one of the initial skyrmion core. As illustrated in the next paragraph, this behaviour can be used for the core-reversal strategy under external applied electric field. The transition from (S) to (P) state goes through some (M) metastable states in a very narrow PMA window.
The variation of the energy density barrier ∆E dens separating the skyrmion core spins and the oppositely oriented edge ones has a non-monotonous evolution with DMI and PMA, as illustrated in Figure 16. In both situations, an optimum window exists for stable skyrmionic states (Figure 16b,c): complex chiral structures with Q > 0 or saturated state (Q = 0) appearing after some up-threshold DMI and PMA values, respectively. A larger PMA implicates a confinement of the skyrmion core in a smaller region. This would require a larger DMI for twisting the spins against the direct exchange interaction. The threshold from (S) to (P) is roughly showing a parabolic increase with the DMI (Figure 15a) dotted white line.  The magnetostatic demagnetizing energy is mainly influenced by two parameters: the saturation magnetization Ms and the size and the shape of the magnetic disks. Both have a direct impact on the DMI-PMA-Q ground state phase diagram, and, implicitly, on the skyrmion stabilization parameters window. This is illustrated by Figure 18   A deeper insight into the skyrmion stability, with respect to the increase in the anisotropy, can be obtained by analyzing the corresponding variation of different contributions to the magnetic free energy ( Figure 17). The analysis has been performed along a path, indicated by a solid black vertical line in Figure 15, corresponding to a fixed value of DMI (D ind = 4 mJ/m 2 ) and variable PMA (K u = 0-2.5 × 10 +6 J/m 3 ). We see that the skyrmion stability is a complex issue, being determined by the balance between the PMA, the DMI, the direct exchange and the demagnetizing (magnetostatic) energies. When the PMA increases during the evolution from a vortex state to a skyrmionic one, the total energy decreases. However, this evolution costs more and more in demagnetizing energy which is continuously increasing. At a certain critical value of PMA, an abrupt transition into a saturated perpendicular state will stop the increase in the magnetostatic energy and will accommodate both the anisotropy and the direct exchange energies, in detriment to the DMI energy (the abrupt transition from (S) to (P) is also illustrated in Figure 17).
The magnetostatic demagnetizing energy is mainly influenced by two parameters: the saturation magnetization M s and the size and the shape of the magnetic disks. Both have a direct impact on the DMI-PMA-Q ground state phase diagram, and, implicitly, on the skyrmion stabilization parameters window. This is illustrated by Figure 18    The magnetostatic demagnetizing energy is mainly influenced by two parameters: the saturation magnetization Ms and the size and the shape of the magnetic disks. Both have a direct impact on the DMI-PMA-Q ground state phase diagram, and, implicitly, on the skyrmion stabilization parameters window. This is illustrated by Figure 18  The influence of the size of the magnetic disks is illustrated in Figure 19, in which we represented calculated DMI-PMA-Q phase diagrams for 90 nm and 50 nm disks patterned from Fe-type materials with M s = 1714 kA/m; A ex = 2.1 × 10 −11 J/m. One can clearly observe that the smaller ferromagnetic disks would need a larger DMI to promote skyrmionic states (because the core should be confined within a smaller area and spins are swirling within a narrow length-scale). A direct consequence would be the necessity to fabricate materials and/or profit from additive effects in complex multilayered architectures to gain access to the large DMI values necessary to host nanometric size skyrmions in nanometric disks, as required for high-density storage or qubit applications. As illustrated in Figure 19c, this issue can be partially overcome by using materials with a larger M s : the minimum value of DMI for skyrmionics ground state stabilization decreases with the increase in M s .
From all the DMI-PMA-Q diagrams (Figures 13a, 15a, 18, 19 and 20) and from the variation of the value of the topological charge Q with the increase in the anisotropy (Figure 17b), one can identify the manifestation of the topological protection of skyrmions. The transition between the skyrmionic domain (|Q| ≈ 1) towards the uniform magnetization one (|Q| ≈ 1) is always abrupt (e.g., when enhancing the anisotropy).
Another investigation that we performed concerned the influence of a perpendicularly applied magnetic field on the magnetic nano-disk ground states. This aspect is illustrated in Figure 20. The first visible/sizeable effect concerns the vanishing of complex chiral states with positive topological charge |Q| > 1 in favor of vortex configurations. Then, the parameter window for stabilizing skyrmions is pushed towards smaller DMI and PMA values and enlarged for the higher PMA range. Moreover, the size of the skyrmion shrinks with increasing the magnetic field (as already shown by Moreau-Luchaire et al. in [9]). This could be an interesting issue for applications from both materials and skyrmion manipulation points of view.
The studies presented in this paragraph relative to the stabilization of different types of chiral ground states illustrate the importance of material engineering, e.g., proper choice and control of M s , A ex , PMA, DMI, to be able to fit within the parameter range required for skyrmion stabilization, in view of any skyrmionic application. These concepts are not only valid for nanostructures constituted from single magnetic layers but can be further extrapolated to more complex heterostructures.  The influence of the size of the magnetic disks is illustrated in Figure 19, in which we represented calculated DMI-PMA-Q phase diagrams for 90 nm and 50 nm disks patterned from Fe-type materials with Ms = 1714 kA/m; Aex = 2.1 × 10 −11 J/m. One can clearly observe that the smaller ferromagnetic disks would need a larger DMI to promote skyrmionic states (because the core should be confined within a smaller area and spins are swirling within a narrow length-scale). A direct consequence would be the necessity to fabricate materials and/or profit from additive effects in complex multilayered architectures to gain access to the large DMI values necessary to host nanometric size skyrmions in nanometric disks, as required for high-density storage or qubit applications. As illustrated in Figure  19c, this issue can be partially overcome by using materials with a larger Ms: the minimum value of DMI for skyrmionics ground state stabilization decreases with the increase in Ms. From all the DMI-PMA-Q diagrams (Figures 13a,15a and 18 -20) and from the variation of the value of the topological charge Q with the increase in the anisotropy (Figure  17b), one can identify the manifestation of the topological protection of skyrmions. The transition between the skyrmionic domain (| | ≈ 1) towards the uniform magnetization one (| | ≈ 1) is always abrupt (e.g., when enhancing the anisotropy).
Another investigation that we performed concerned the influence of a perpendicularly applied magnetic field on the magnetic nano-disk ground states. This aspect is illus-  The studies presented in this paragraph relative to the stabilization of different types of chiral ground states illustrate the importance of material engineering, e.g., proper choice and control of Ms, Aex, PMA, DMI, to be able to fit within the parameter range required for skyrmion stabilization, in view of any skyrmionic application. These concepts are not only valid for nanostructures constituted from single magnetic layers but can be further extrapolated to more complex heterostructures.
Indeed, in Figure 21, we present an example of a DMI-PMA-Q ground state phase diagram corresponding to an antiferromagnetic (AF) configuration. The diagram looks very similar to the one corresponding to a single layer, the chiral structures (e.g., skyrmions) in the two layers are antiferromagnetically mirrored and have opposite signs of the topological charge (Figure 21c). Such kinds of skyrmions could be stabilized in synthetic antiferromagnetic multilayers. They would have the advantage that the two FM layers which can be finely tuned are joined together for an increased thermal stability and improved dynamical behavior of the Skyrmion, e.g., correction of the "Magnus" drift or skyrmionic spin Hall effect (Figure 22), observed if single skyrmions move in racetracks [31]. Moreover, as recently suggested, such kind of AF-coupled skyrmions could be used as coupled qubits [19]. Indeed, in Figure 21, we present an example of a DMI-PMA-Q ground state phase diagram corresponding to an antiferromagnetic (AF) configuration. The diagram looks very similar to the one corresponding to a single layer, the chiral structures (e.g., skyrmions) in the two layers are antiferromagnetically mirrored and have opposite signs of the topological charge ( Figure 21c). Such kinds of skyrmions could be stabilized in synthetic antiferromagnetic multilayers. They would have the advantage that the two FM layers which can be finely tuned are joined together for an increased thermal stability and improved dynamical behavior of the Skyrmion, e.g., correction of the "Magnus" drift or skyrmionic spin Hall effect ( Figure 22), observed if single skyrmions move in racetracks [31]. Moreover, as recently suggested, such kind of AF-coupled skyrmions could be used as coupled qubits [19].
At the end of this subsection, we would like to illustrate that multi-skyrmionic ground states can be also generated in nano-disks by properly tuning the M s , PMA, and DMI parameters ( Figure 23). The simulations presented in this figure are performed considering a 250 nm diameter disk, large enough to be able to accommodate multiple skyrmions within a PMA-DMI window in the range of phase diagram parameters of Figure 14c; the simulations being for a CoFeB type material with M s = 580 kA/m; A ex = 1.5 × 10 −11 J/m.
The possibility to generate such kind of multiple skyrmionic states by PMA and DMI tuning (e.g., by electric field gating), is particularly important for memristor functions, e.g., in a racetrack, where the resistance read by a magnetic tunnel junction depends on the number of skyrmions displaced by the current [18].

Skyrmionic Nano-Oscillators
Once the skyrmions are generated, an important issue for a different type of applications is their manipulation in a different type of device.  At the end of this subsection, we would like to illustrate that multi-skyrmionic ground states can be also generated in nano-disks by properly tuning the Ms, PMA, and DMI parameters (Figure 23). The simulations presented in this figure are performed considering a 250 nm diameter disk, large enough to be able to accommodate multiple skyrmions within a PMA-DMI window in the range of phase diagram parameters of Figure  14c; the simulations being for a CoFeB type material with Ms = 580 kA/m; Aex = 1.5 × 10 −11 J/m.  At the end of this subsection, we would like to illustrate that multi-skyrmionic ground states can be also generated in nano-disks by properly tuning the Ms, PMA, and DMI parameters (Figure 23). The simulations presented in this figure are performed considering a 250 nm diameter disk, large enough to be able to accommodate multiple skyrmions within a PMA-DMI window in the range of phase diagram parameters of Figure  14c; the simulations being for a CoFeB type material with Ms = 580 kA/m; Aex = 1.5 × 10 −11 J/m. The T pulse -J c -Q and T pulse -J c -m z phase diagrams illustrate that for writing a skyrmion by STT, a spin-polarized current pulse must be applied, with properly chosen J c and T pulse . After cutting the pulse, the skyrmionic nucleated pattern will have a damped oscillation whose size and shape is evanescently converging towards the final diameter imposed by the magnetic material parameters (DMI, PMA, M s ) and the patterned disk geometry. However, as if for a mechanical damped oscillator, we would supply energy to compensate the damping losses, e.g., by acting with an external driving torque provided by a spin polarized current that competes with the damping torque, the skyrmion can be driven in a steady self-sustained oscillation regime: in which the size of the skyrmion is cyclically increasing and decreasing (so-called "breathing" regime). In Figure 24, we illustrate micromagnetically simulated magnetization dynamics feature for the writing (damped oscillations from pulse nucleated pattern towards the final skyrmionic state) and the steady oscillation regime. They concern a skyrmion generated in a 90 nm diameter nano-disk, patterned from an Au/Fe(5ML)/MgO multilayer system with the magnetic parameters extracted from ab-initio calculations: K u (PMA) = 1.91 × 10 +6 J/m 3  The possibility to generate such kind of multiple skyrmionic states by PMA and DMI tuning (e.g., by electric field gating), is particularly important for memristor functions, e.g., in a racetrack, where the resistance read by a magnetic tunnel junction depends on the number of skyrmions displaced by the current [18].

Skyrmionic Nano-Oscillators
Once the skyrmions are generated, an important issue for a different type of applications is their manipulation in a different type of device.
The Tpulse-Jc-Q and Tpulse-Jc-mz phase diagrams illustrate that for writing a skyrmion by STT, a spin-polarized current pulse must be applied, with properly chosen Jc and Tpulse. After cutting the pulse, the skyrmionic nucleated pattern will have a damped oscillation whose size and shape is evanescently converging towards the final diameter imposed by the magnetic material parameters (DMI, PMA, Ms) and the patterned disk geometry. However, as if for a mechanical damped oscillator, we would supply energy to compensate the damping losses, e.g., by acting with an external driving torque provided by a spin polarized current that competes with the damping torque, the skyrmion can be driven in a steady self-sustained oscillation regime: in which the size of the skyrmion is cyclically increasing and decreasing (so-called "breathing" regime). In Figure 24, we illustrate micromagnetically simulated magnetization dynamics feature for the writing (damped oscillations from pulse nucleated pattern towards the final skyrmionic state) and the steady oscillation regime. They concern a skyrmion generated in a 90 nm diameter nano-disk, patterned from an Au/Fe(5ML)/MgO multilayer system with the magnetic parameters extracted from ab-initio calculations: Ku (PMA) = 1.91 × 10 +6 J/m 3   The conditions used for writing and sustaining the oscillations were: Jc = 1 × 10 +12 A/m 2 and Tpulse = 1.5 ns for writing from an initial P state of the MTJ nanopillar, then Jc = −1.5 × 10 +12 A/m 2 for STT driving torque driven oscillations. The spin-polarized current Jc injected into the nanopillar has a double effect: the damping-like component will provide the negative torque for sustaining the steady oscillations and the field-like component will modulate the frequency of oscillation determining a strongly nonlinear variation with the current density Jc. For maintaining the steady oscillation regime of the skyrmionic nanooscillator, Jc must be adjusted in a proper range. If it is too small, the damping will not be The conditions used for writing and sustaining the oscillations were: J c = 1 × 10 +12 A/m 2 and T pulse = 1.5 ns for writing from an initial P state of the MTJ nanopillar, then J c = −1.5 × 10 +12 A/m 2 for STT driving torque driven oscillations. The spin-polarized current J c injected into the nanopillar has a double effect: the damping-like component will provide the negative torque for sustaining the steady oscillations and the field-like component will modulate the frequency of oscillation determining a strongly nonlinear variation with the current density J c . For maintaining the steady oscillation regime of the skyrmionic nano-oscillator, J c must be adjusted in a proper range. If it is too small, the damping will not be entirely compensated, and the oscillations will gradually vanish. J c values which are too large will provide anti-damping torques which are also too large. The skyrmion core will increase up to the limit of the disk size, and collapse in a stable saturated state. The skyrmionic nano-oscillators [32] gained a lot of interest in recent years. Their wide window of operating frequency, in the range of the microwaves (1-100 GHz), tunable by external stimuli (magnetic fields, STT of spin-polarized currents), open interesting application perspectives as nanoscale electric oscillators, and sensitive magnetic field sensors [33]. Their strongly nonlinear dynamics can be used in studies of chaotic phenomena [34] and for neuro-inspired devices suitable for neuromorphic applications, e.g., pattern recognition [35].

Skyrmion Manipulation by Electric Fields
Once the skyrmion is generated, the next issue is to manipulate them in skillfully designed classic (storage, logics), neuromorphic and quantum devices. The strategy of skyrmion manipulation by electric fields is one of the most energetically efficient options. It is based on the electric field control of the perpendicular magnetic anisotropy [36]. In Figure 25a, we zoom on a skyrmionic zone within a ground states DMI-PMA-Q diagram. We set an initial skyrmionic state (i), described by PMA = 1.675 × 10 +6 J/m 3 and DMI = 2.4 × 10 −3 J/m 2 , chosen to be in the proximity of the transition zone between the skyrmionic (S) and the saturated (P) states. This initial skyrmion, being nucleated from a saturated state along −Oz m = (0,0,−1), has a negative core polarization, the other parameters used in the simulation are M s = 1714 × 10 +3 A/m, A ex = 2.1 × 10 −11 J/m and α = 0.01 (Fe). Afterwards, we simulate the application of an electric field pulse with a properly adjusted delay (here 0.3 ns) and with an orientation providing an enhancement of the PMA. The E-field increase in PMA would lead to the progressive decrease in the skyrmion diameter, up to the complete erasure (Figure 25b states (1)-(5)). After cutting down the pulse, the anisotropy is restored to its initial value and a new skyrmionic state with an Afterwards, we simulate the application of an electric field pulse with a properly adjusted delay (here 0.3 ns) and with an orientation providing an enhancement of the PMA. The E-field increase in PMA would lead to the progressive decrease in the skyrmion diameter, up to the complete erasure (Figure 25b states (1)-(5)). After cutting down the pulse, the anisotropy is restored to its initial value and a new skyrmionic state with an opposite core orientation will be relaxed as the ground state. In Figure 25c, we illustrate the line profiles of the total magnetic energy density corresponding to the states depicted in Figure 25a. The core polarization always corresponds to the orientation of the disk magnetization in the saturated state before relaxation. When erasing a skyrmion with negative polarization by core shrinkage, we end with an up-saturated state and then, the new relaxed skyrmionic ground state will have a positive core polarization, and vice-versa.
Following this strategy, as illustrated in Figure 26, the skyrmionic core can be cyclically reversed from up to down orientation. These results demonstrate the capability of the electric field control of a binary (0)/(1) information that could be encoded in the skyrmion core magnetization, e.g., if the magnetic disk containing the skyrmion is the free magnetic layer of a magnetic tunnel junction whose resistance depends on the skyrmionic core polarization. On the other hand, an interesting issue would be related to the possibility to switch the chirality and the core polarization of a skyrmion, simultaneously. Our simulations indicate that this can be performed by exploiting the electric field modulation of the DMI. From ab-initio calculations, we can determine the rate of the electric field variation of the DMI, as illustrated in Figure 27a. We fix an initial state (i) for a value of DMI allowing a ground state skyrmionic stabilization (blue open circle). Then, from the theoretical DMI vs E-field curve, the value (and the orientation) of the electric field required to reverse and obtain a positive DMI value, providing a new (f) skyrmionic ground state with opposite chirality, can be extrapolated. Following the path between (i) and (f) the simulated magnetization dynamics show that the initial and the final skyrmion have opposite core polarization and chirality (as in the case of antiferromagnetically mirroring). This possibility for core and chirality reversal can be useful in storage or logic devices where, by E-field On the other hand, an interesting issue would be related to the possibility to switch the chirality and the core polarization of a skyrmion, simultaneously. Our simulations indicate that this can be performed by exploiting the electric field modulation of the DMI. From ab-initio calculations, we can determine the rate of the electric field variation of the DMI, as illustrated in Figure 27a. We fix an initial state (i) for a value of DMI allowing a ground state skyrmionic stabilization (blue open circle). Then, from the theoretical DMI vs E-field curve, the value (and the orientation) of the electric field required to reverse and obtain a positive DMI value, providing a new (f) skyrmionic ground state with opposite chirality, can be extrapolated. Following the path between (i) and (f) the simulated magnetization dynamics show that the initial and the final skyrmion have opposite core polarization and chirality (as in the case of antiferromagnetically mirroring). This possibility for core and chirality reversal can be useful in storage or logic devices where, by E-field gating one can change the sign of the skyrmionic Hall effect and the direction of the skyrmionic drift. However, as illustrated in Figure 27a, if only intrinsic effects are responsible for the DMI variation with the electric field, a quite large gating voltage is required. From a practical point of view, this could bring limitations due to the dielectric breakdown issues of the insulator used for the voltage gating (experimental geometry depicted in the top insert of Figure 27). Therefore, beyond optimizing the quality of the dielectrics used for the voltage gating, a major challenge for the implementation of the DMI sign reversal in skyrmionic applications would be to design and elaborate materials and related multilayered configurations in which the skyrmionic states can be stabilized at low DMI and PMA values.

Multiscale Modeling of Skyrmionic Nanomaterials
The analysis illustrated in Figure 27 demonstrates that a multiscale strategy, combining ab-initio and micromagnetic modeling, can be particularly useful for theoretically investigating skyrmionic or other topological spin texture materials. In our case, the ab-initio analysis provided a fundamental insight into the PMA and DMI mechanisms and quantified their electric field control capability in various ab-initio materials and multilayer architectures experimentally conceivable by standard deposition techniques (i.e., Molecular Beam Epitaxy, Sputtering). As sketched in Figure 28, this would indicate the theoretically possible paths of "motion" within a DMI-PMA-Q diagram to access and manipulate the available different chiral magnetic states. Extrapolated beyond this example, in the first stage, the ab-initio analysis could be employed to get insight into the anatomy of intrinsic properties and phenomena allowing to theoretically design materials or multilayer configurations with optimum magnetic properties. Then, in the next stage, their mesoscopic static and dynamic magnetic properties can be further tested/simulated using micromagnetic tools. The multiscale modeling can be also done the other way around. From micromagnetic extended phase diagrams one can identify the range of optimum However, as illustrated in Figure 27a, if only intrinsic effects are responsible for the DMI variation with the electric field, a quite large gating voltage is required. From a practical point of view, this could bring limitations due to the dielectric breakdown issues of the insulator used for the voltage gating (experimental geometry depicted in the top insert of Figure 27). Therefore, beyond optimizing the quality of the dielectrics used for the voltage gating, a major challenge for the implementation of the DMI sign reversal in skyrmionic applications would be to design and elaborate materials and related multilayered configurations in which the skyrmionic states can be stabilized at low DMI and PMA values.

Multiscale Modeling of Skyrmionic Nanomaterials
The analysis illustrated in Figure 27 demonstrates that a multiscale strategy, combining ab-initio and micromagnetic modeling, can be particularly useful for theoretically investigating skyrmionic or other topological spin texture materials. In our case, the abinitio analysis provided a fundamental insight into the PMA and DMI mechanisms and quantified their electric field control capability in various ab-initio materials and multilayer architectures experimentally conceivable by standard deposition techniques (i.e., Molecular Beam Epitaxy, Sputtering). As sketched in Figure 28, this would indicate the theoretically possible paths of "motion" within a DMI-PMA-Q diagram to access and manipulate the available different chiral magnetic states. Extrapolated beyond this example, in the first stage, the ab-initio analysis could be employed to get insight into the anatomy of intrinsic properties and phenomena allowing to theoretically design materials or multilayer configu-rations with optimum magnetic properties. Then, in the next stage, their mesoscopic static and dynamic magnetic properties can be further tested/simulated using micromagnetic tools. The multiscale modeling can be also done the other way around. From micromagnetic extended phase diagrams one can identify the range of optimum parameters for the envisaged magnetic configuration and nanoscopic property. Then, ab-initio calculations can be used to design the magnetic material or complex nanostructure that would provide those properties. If such a kind of multiscale analysis precedes the experimental elaboration of complex nanomaterials, the net advantage would be the reduction in the experimental efforts, costs, and time consumption. that would provide those properties. If such a kind of multiscale analysis precedes the experimental elaboration of complex nanomaterials, the net advantage would be the reduction in the experimental efforts, costs, and time consumption.

Experimental Issues on Magnetic Skyrmionic Nano-Materials
Based on the predictive theoretical calculations, we got insight into the magnetic properties of selected materials and layered architectures and placed them on the different types of phase diagrams describing different magnetic chiral states. Then, we illustrated the possibility to manipulate the magnetic skyrmionic states for various applications. The last issue that we would like to discuss is about some currently available magnetic materials and corresponding layered thin film architectures that could accommodate the DMI, PMA and Ms parameter range of our calculated phase diagrams. A selection with recent data available from the literature, concerning selected values of the saturation magnetization, DMI and PMA is illustrated by Figure 29 (the corresponding references are indicated). Concerning the saturation magnetization of a ferromagnetic layer - Figure 29bconsidering its ultrathin thickness range, the surface contributions will be significant. Therefore, the value of Ms depends on the thin film thickness and the chemical nature of the top (X) and bottom (Y) interface layers X/FM/Y. From Figure 29a, one can see that the anisotropy range is well covered by available experimental systems. They not only provide a wide range of surface anisotropies, but one also has the possibility to tune the value of Ku by the thickness of the ferromagnetic layer: Ku=Ks/t. On the other hand, the currently available DMI range, for the most part of the studied systems, is beyond 3 mJ/m 2 . As mentioned, the DMI can be enhanced by additive effects [9] or skillful engineering of multilayered sequences, use of synthetic antiferromagnets [37], etc. Indeed, recently, in the epitaxial W(110)/Fe/Co bilayer system a giant DMI, of 6.55 J/m 2 obtained by quantum engineering of the lattice symmetry has been demonstrated [38]. On the other hand, in this paper we show that ab-initio and micromagnetic complementary tools can be successfully used to calculate the magnetic properties of multilayer skyrmionic materials, the predicted DMI and PMA values being similar with those measured in experimental stacks. Moreover, we theoretically indicate a path for enhancing the DMI, the PMA and their response to external electric fields by the skillful engineering of interfaces: a monolayer of

Experimental Issues on Magnetic Skyrmionic Nano-Materials
Based on the predictive theoretical calculations, we got insight into the magnetic properties of selected materials and layered architectures and placed them on the different types of phase diagrams describing different magnetic chiral states. Then, we illustrated the possibility to manipulate the magnetic skyrmionic states for various applications. The last issue that we would like to discuss is about some currently available magnetic materials and corresponding layered thin film architectures that could accommodate the DMI, PMA and M s parameter range of our calculated phase diagrams. A selection with recent data available from the literature, concerning selected values of the saturation magnetization, DMI and PMA is illustrated by Figure 29 (the corresponding references are indicated). Concerning the saturation magnetization of a ferromagnetic layer - Figure 29b-considering its ultrathin thickness range, the surface contributions will be significant. Therefore, the value of M s depends on the thin film thickness and the chemical nature of the top (X) and bottom (Y) interface layers X/FM/Y. From Figure 29a, one can see that the anisotropy range is well covered by available experimental systems. They not only provide a wide range of surface anisotropies, but one also has the possibility to tune the value of K u by the thickness of the ferromagnetic layer: K u = K s /t. On the other hand, the currently available DMI range, for the most part of the studied systems, is beyond 3 mJ/m 2 . As mentioned, the DMI can be enhanced by additive effects [9] or skillful engineering of multilayered sequences, use of synthetic antiferromagnets [37], etc. Indeed, recently, in the epitaxial W(110)/Fe/Co bilayer system a giant DMI, of 6.55 J/m 2 obtained by quantum engineering of the lattice symmetry has been demonstrated [38]. On the other hand, in this paper we show that ab-initio and micromagnetic complementary tools can be successfully used to calculate the magnetic properties of multilayer skyrmionic materials, the predicted DMI and PMA values being similar with those measured in experimental stacks. Moreover, we theoretically indicate a path for enhancing the DMI, the PMA and their response to external electric fields by the skillful engineering of interfaces: a monolayer of Pt inserted at the top Fe/MgO would enhance the PMA by almost an order of magnitude, from 1.91 to 15.77 MJ/m 3 , the DMI from 1.67 to 2.27 mJ/m 2 and their response to an external electric field roughly by a factor of 2.
However, the yet existing experimental "gap", opening roughly above 3 mJ/m 2 , shows that there is still plenty of space for theoretical and experimental research concerning the topological spin texture materials. A promising issue would be to combine interfacial and bulk DMI mechanisms requiring magnetic materials with intrinsic large spin-orbit interactions. Rare-earth transition metal alloys are more and more considered promising candidates for small skyrmions and ultrafast chiral spin texture dynamics applications [39]. Within this emerging topic, skyrmionic bubbles have been recently obtained [40] in Rare Earth (RE)-based REMn 2 Ge 2 (RE = Ce, Pr, Nd) magnets, in a wide temperature range (220-320 K). Lattices of magnetic skyrmions have been observed in centrosymmetric rare earth compounds, such as Gd 2 PdSi 3 and GdRu 2 Si 2 [41]. Spontaneous topological magnetic transitions were identified in NdCo 5 RE magnets [42]. Compact ferrimagnetic skyrmions were observed in DyCo 3 films [43] and interfacial chiral magnetism and isolated skyrmions have been demonstrated in SmCo 5 -based magnetic multilayers featuring perpendicular magnetic anisotropy [44].
The topological chiral textures' magnetic material spectrum is not only restricted to metals but also opened to insulating magnetic oxides. Recently, it was demonstrated that the DMI can be significantly controlled by rare earth orbital magnetism and strain engineering in insulating magnetic oxides [45]. Even the paradigm that skyrmion occurrence is strictly related to the lack of symmetry centers in crystals or interfaces has been newly invalidated by showing that skyrmions can be created in materials having centers of symmetry but possessing geometric frustration [46]. However, the yet existing experimental "gap", opening roughly above 3 mJ/m 2 , shows that there is still plenty of space for theoretical and experimental research concerning the topological spin texture materials. A promising issue would be to combine interfacial and bulk DMI mechanisms requiring magnetic materials with intrinsic large spinorbit interactions. Rare-earth transition metal alloys are more and more considered promising candidates for small skyrmions and ultrafast chiral spin texture dynamics applications [39]. Within this emerging topic, skyrmionic bubbles have been recently obtained [40] in Rare Earth (RE)-based REMn2Ge2 (RE = Ce, Pr, Nd) magnets, in a wide temperature range (220-320 K). Lattices of magnetic skyrmions have been observed in centrosymmetric rare earth compounds, such as Gd2PdSi3 and GdRu2Si2 [41]. Spontaneous topological magnetic transitions were identified in NdCo5 RE magnets [42]. Compact ferrimagnetic skyrmions were observed in DyCo3 films [43] and interfacial chiral magnetism and isolated skyrmions have been demonstrated in SmCo5-based magnetic multilayers featuring perpendicular magnetic anisotropy [44].
The topological chiral textures' magnetic material spectrum is not only restricted to metals but also opened to insulating magnetic oxides. Recently, it was demonstrated that the DMI can be significantly controlled by rare earth orbital magnetism and strain engineering in insulating magnetic oxides [45]. Even the paradigm that skyrmion occurrence is strictly related to the lack of symmetry centers in crystals or interfaces has been newly invalidated by showing that skyrmions can be created in materials having centers of symmetry but possessing geometric frustration [46]. When only the surface contribution to the PMA (Ks) was available, the anisotropy Ku has been calculated considering a ferromagnetic thickness layer of t = 1 nm: Ku = Ks/t; when the anisotropy field = (T) was available for a layer with given Ms (A/m), Ku (J/m 3 ) was calculated as = /2 [47,48].
Within these complex and challenging problematics, multiscale studies combining theoretical and experimental investigations could complete the missing elements in the puzzle and answer the complex applicative issues.

Conclusions
Within a multiple-scale modeling framework, we addressed the problems of magnetic skyrmions. On specially designed multilayer stacks, ab-initio calculations provided insight into characteristic magnetic quantities of multi-layered heterostructures: the saturation magnetization, the anisotropy, the DMI and their response to a gating electric field. When only the surface contribution to the PMA (K s ) was available, the anisotropy K u has been calculated considering a ferromagnetic thickness layer of t = 1 nm: K u = K s /t; when the anisotropy field B k = µ o H k (T) was available for a layer with given M s (A/m), K u (J/m 3 ) was calculated as K u = M s H k /2 [47,48].
Within these complex and challenging problematics, multiscale studies combining theoretical and experimental investigations could complete the missing elements in the puzzle and answer the complex applicative issues.

Conclusions
Within a multiple-scale modeling framework, we addressed the problems of magnetic skyrmions. On specially designed multilayer stacks, ab-initio calculations provided insight into characteristic magnetic quantities of multi-layered heterostructures: the saturation magnetization, the anisotropy, the DMI and their response to a gating electric field. Using as input these results, micromagnetic calculations solving the Landau-Lifshitz-Gilbert dynamics, allowed us to identify critical regimes for writing skyrmions in patterned nanopillars. Extended phase diagrams, in which the magnetic textures were classified according to their topological charge, allowed us to identify the suitable range of magnetic anisotropy, DMI and saturation magnetization that allow the stabilization of skyrmionic ground states. An analysis of different contributions to the total magnetic free energy underlined critical issues influencing the skyrmions' stability. Then, some manipulation strategies of the skyrmions' core polarization and chirality by electric fields have been illustrated. Finally, we located the simulation window of parameters used in our extended modeling within the current experimental reality, in terms of existing or the perspective magnetic materials compatible with skyrmionic applications in classic, neuromorphic and quantum information technologies.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A
The Influence of Thermal Fluctuations at Finite Temperature The ab-initio as well as the micromagnetic calculations presented so far have been performed at zero Kelvin. This simplified approach, which is not considering the thermal blur effects, allowed us to clearly observe and critically predict critical the magnetic properties and behavior. However, for the experimental extrapolation of the results, one may wonder about the thermal stability of the theoretically predicted chiral structures (e.g., skyrmions). To answer to this question, we have also performed sets of micromagnetic simulations at finite temperature. In the Mumax3 code, the temperature is implemented by means of a fluctuating thermal field: → B therm = → η (step) 2µ 0 αk B T sim γ 0 µ 0 M s V sim ∆t , where α -damping parameter, k B -Boltzmann constant, T sim -simulation temperature, γ 0 -gyromagnetic ratio, µ 0 -magnetic permeability, M s -saturation magnetization, V sim -the simulation volume (cell) and ∆tthe simulation time step, → η (step) is a random vector from a standard normal distribution whose value is changed every step.
As illustrated in Figure A1, with increasing the temperature, the thermal field is randomly disordering the local spins within a certain chiral structure. Figure A1, panel (a) shows the DMI-PMA-Q phase diagrams calculated at 0 K, 100 K and 300 K. The panel (b) illustrates a skyrmionic structure, corresponding to the white point in each diagram, for each temperature. At a finite temperature, locally, the spin orientation will have a random distribution around the "equilibrium" zero-Kelvin position. The larger the temperature, the larger the width of this distribution. The temperature and DMI dependent random distribution of local spins orientation will further affect the value of the numerically calculated topological charge: = ⃗ • ⃗ × ⃗ , as illustrated in Figure A1, panel (a).
A qualitative illustration concerning the stability in temperature of a skyrmionic structure, is also presented in Figure A2, in which vector field glyph representation is superimposed to corresponding scalar field mz maps, for each temperature. At a finite temperature, locally, the spin orientation will have a random distribution around the "equilibrium" zero-Kelvin position. The larger the temperature, the larger the width of this distribution. The temperature and DMI dependent random distribution of local spins orientation will further affect the value of the numerically calculated topological dy , as illustrated in Figure A1, panel (a). A qualitative illustration concerning the stability in temperature of a skyrmionic structure, is also presented in Figure A2, in which vector field glyph representation is superimposed to corresponding scalar field m z maps, for each temperature.
One can observe that, for the chosen value of the DMI and PMA, the skyrmionic structure survives even at room temperature (T = 300 K). However, for temperatures above 300 K, the skyrmionic state will finally be erased by the local thermal field. Because the barrier protecting the spins from the core to reverse their orientation is controlled by the DMI (as illustrated in Figure 12b), the thermal stability of the skyrmions is increasing with the DMI increase (as qualitatively illustrated in Figure A3), or by applying a magnetic field.
local spins orientation will further affect the value of the numerically calculated topological charge: = ⃗ • ⃗ × ⃗ , as illustrated in Figure A1, panel (a).
A qualitative illustration concerning the stability in temperature of a skyrmionic structure, is also presented in Figure A2, in which vector field glyph representation is superimposed to corresponding scalar field mz maps, for each temperature. One can observe that, for the chosen value of the DMI and PMA, the skyrmionic structure survives even at room temperature (T = 300 K). However, for temperatures above 300K, the skyrmionic state will finally be erased by the local thermal field. Because the barrier protecting the spins from the core to reverse their orientation is controlled by the DMI (as illustrated in Figure 12b), the thermal stability of the skyrmions is increasing with the DMI increase (as qualitatively illustrated in Figure A3), or by applying a magnetic field. Figure A3. Skyrmionic state at T = 300 K for: (a) DMI = 5 mJ/m 2 and (b) DMI = 2.5 mJ/m 2 . The calculation has been performed for a 90 nm disk with the following parameters: DMI = 5 mJ/m 2 , Ku = 1.7 × 10 +6 J/m 3 , Ms = 1714 kA/m, Aex = 2.1 × 10 +11 J/m; the Gilbert damping constant α = 0.01. In each panel, the top image is a glyph representation of the magnetization vector field within the disk. The bottom image is a scalar field representation of the mz values.