Encoding and Multiplexing Information Signals in Magnetic Multilayers with Fractional Skyrmion Tubes

Tailored magnetic multilayers (MMLs) provide skyrmions with enhanced thermal stability, leading to the possibility of skyrmion-based devices for room-temperature applications. At the same time, the search for additional stable topological spin textures has been under intense research focus. Besides their fundamental importance, such textures may expand the information encoding capability of spintronic devices. However, fractional spin texture states within MMLs in the vertical dimension are yet to be investigated. In this work, we demonstrate numerically fractional skyrmion tubes (FSTs) in a tailored MML system. We subsequently propose to encode sequences of information signals with FSTs as information bits in a tailored MML device. Micromagnetic simulations and theoretical calculations are used to verify the feasibility of hosting distinct FST states within a single device, and their thermal stability is investigated. A multilayer multiplexing device is proposed, where multiple sequences of the information signals can be encoded and transmitted based on the nucleation and propagation of packets of FSTs. Finally, pipelined information transmission and automatic demultiplexing are demonstrated by exploiting the skyrmion Hall effect and introducing voltage-controlled synchronizers and width-based track selectors. The findings indicate that FSTs can be potential candidates as information carriers for future spintronic applications.


INTRODUCTION
Magnetic skyrmions are particle-like topological spin configurations. 1,2 They are promising candidates as information carriers in future information technologies, owning to their topological protection arising from their unique spin texture and dynamical properties. 2,3 This stability comes from the competition between magnetic interactions which favor colinear spin configurations, such as Heisenberg exchange coupling and perpendicular magnetic anisotropy (PMA), and interactions favoring orthogonal configurations such as dipolar coupling and Dzyaloshinskii−Moriya interaction (DMI). 4,5 In both bulk magnetic and magnetic multilayer (MML) systems with broken inversion symmetry, skyrmion spin states become one of the system energy minima. The magnetic skyrmion can be described by an integer topological index, called the skyrmion number N, which counts how many times the magnetization wraps around a unit sphere. The skyrmion number is defined as 6 (1) where m is the normalized magnetization and N = ±1 corresponds to the case of magnetic skyrmions, with the sign reflecting the polarity. Skyrmions can be competitive candidates as information carriers in low power and highly efficient computational devices because of their non-volatility, nanoscale size, and ease of manipulation. 3 These advantages have inspired proposals for their implementation in skyrmionic transistors, 7 logic gates, 8−10 racetrack memory, 11−13 nanooscillators, 14 resonant diodes, 15 neuromorphic computing, 16,17 and reservoir computing. 18 Notably, ref 17 is an experimental demonstration of a skyrmionic synapse device, while the other references contain numerical results.
For each of these applications, room-temperature operation is a critical requirement for realistic device integration. Recently, tailored MMLs have been explored as a means to host skyrmions at room temperature, 19,20 where the stacking of a repeated layer structure leads to enhanced thermal stability. This stability in MMLs can be attributed to the increased DMI from the asymmetric interfaces [heavy metal (HM 1 )| ferromagnetic (FM)] and [FM|HM 2 ] and the increased magnetic volume of skyrmions. Additional to the requirement for room-temperature stability, in order for skyrmions to be practicably used in functional devices it is highly desirable for them to be stabilized under zero applied magnetic field. Recent work has highlighted this possibility, for example, by leveraging lateral confinement effects 21 or by careful engineering of interlayer exchange coupling. 22 At the same time, there has been much recent effort to find skyrmion-like quasiparticles with topological charge other than ±1. Such skyrmionic quasiparticles have been proposed as information carriers to enhance device functionality. 23,24 For example, skyrmionium, which comprises a skyrmion−antiskyrmion pair, is a type of magnetic quasiparticle with a vanishing topological charge N = 0 which has been proposed for use in racetrack memory applications. 23 Similarly, antiskyrmionite (i.e., a double-antiskyrmion−skyrmion pair with N = ∓1 25 ) has the potential to be used as an additional information carrier in skyrmionic devices. More recently, several studies have shown that nanomagnets can indeed host a plethora of topological and non-topological quasiparticles, including both theoretical calculations 24 and experimental demonstrations of skyrmion bags in liquid crystals 26 and skyrmion bundles in chiral magnets. 27 However, these proposals are explorations of skyrmionic particles limited in two-dimensional systems. A natural next step is to explore skyrmion-like quasiparticles which are distinguishable by their three-dimensional profile, including variations in the filmperpendicular direction. In fact, partially stabilized skyrmion quasiparticles have recently been theoretically predicted in bulk chiral magnets, including chiral bobbers 28 and dipole strings, 29 enriching the diversity of the skyrmionic family.
Previous studies in MML systems have mainly focused on skyrmion "tube" states, in which skyrmions are stabilized throughout the MML stack. Standard tubular skyrmions are effective stacks of skyrmions in neighboring magnetic layers which exist across the whole multilayer structure. 20 There has been much recent interest in exploring topological magnetic textures in three dimensions by engineering the properties of heterostructures so as to host complex depth-dependent states. 30 Our proposal is to use fractional skyrmion tubes (FSTs) in MMLs to encode multiple information signals. These FSTs comprise coupled skyrmions that occupy discrete fractions of the full multilayer system. This approach enables us to investigate skyrmion tubes with distinct depth profile, Magnetic parameters in simulations of (c,d) external magnetic field 90 mT and DMI constant 1.9 mJ m −2 which are selected from (b) in order to obtain stable FST states in the system. making it possible to explore multiple states within a single system. We will show that such FSTs can exhibit effective tunability in their thermal stability. Based on their topological properties and propagation by electric currents, we propose a multilayer multiplexing device which relies on the nucleation, propagation, and automatic selection of multiple distinct FSTs in a single device. Our proposal in this work highlights the potential of encoding information via distinct FST states in MMLs. Specifically, we use simulations to demonstrate pipelined information transmission and automatic demultiplexing of information signals via use of voltage-controlled synchronizers and track selectors, respectively. Most results of this work are performed via the micromagnetic simulation package Mumax 3 . 31 Mathematical calculations are performed via the Python mathematical module SciPy. 32

FSTs in MMLs as Multi-Bit Information Carriers.
2.1.1. Potential Multi-Bit Information Carriers. As aforementioned, tailored MMLs have been proposed to stabilize room-temperature skyrmions, typically as skyrmion tubes which extend vertically throughout the depth of the MML. 33,34 In contrast, here we investigate the possibility of obtaining fractions of full skyrmion tubes, wherein skyrmions are stabilized in select layers of the system. All FST states studied in this work comprise Neél-type skyrmions in individual layers, which are favored over Bloch-type skyrmions in multilayer structures with interfacial DMI. Figure 1a shows an example of a 4MML nanotrack hosting four distinct skyrmion states, each comprising a stack of skyrmions which extend over one, two, three, or all four layers of the system. Such skyrmion states are hereafter referred to as FSTs. Thus, Figure 1a exhibits four FSTs, namely (from left to right) a 4MML FST, a 1MML FST, a 3MML FST, and a 2MML FST. We will show that each of these can be nucleated and manipulated individually, leading to a fourfold enhancement of the capacity of the 4MML device to store and process information.
To demonstrate the effect of the micromagnetic parameters on the stability of the four distinct FST states in our example 4MML system, we conducted micromagnetic simulations which explore the phase diagram of FSTs by modifying (i) the external out-of-plane magnetic field, (ii) the DMI constant, and (iii) the FM interlayer exchange coupling. The interlayer exchange coupling modeled throughout this work is of Ruderman−Kittel−Kasuya−Yosida (RKKY) type and is always positive so as to favor parallel alignment between spins in neighboring layers. In physical multilayered systems, the magnitude and sign of the interlayer exchange coupling has an oscillatory dependence on the magnetic layer separation, with the exact character of this sensitivity depending on the particular material structure. Therefore, the interlayer exchange coupling can be effectively modified by tuning the thickness or changing the spacer layer material, 35−37 where either FM or antiferromagnetic interlayer coupling can be achieved. In this work, we model changes to the interlayer exchange coupling without consideration for any particular material system and therefore do not model the corresponding changes to material thickness. To model the RKKY-type interlayer exchange interaction in the micromagnetic simulations, we implement a custom field into the Mumax 3 code which couples nearest neighboring magnetic layers. Next-nearest neighbor and higher order interactions between the magnetic layers are assumed to be negligible. The phase diagrams of distinct FST states in a 4MML film with J Interlayer = 0 and J Interlayer = 0.04 mJ m −2 are shown in Figure 1b.
In the simulations, the external magnetic field was scanned in the range of 10−100 mT with a 10 mT step applied out-ofplane, and the DMI constants were varied from 1.5 to 2.6 mJ m −2 with a step of 0.1 mJ m −2 , while keeping the interlayer exchange coupling constant at 0 and 0.04 mJ m −2 . For each data point in Figure 1b, we configured the system with an initial ansatz that contained a 1MML FST, a 2MML FST, a 3MML FST, and a 4MML FST, respectively. We then allowed the system to equilibrate with the corresponding initial states and magnetic parameters. Each cell of Figure 1b shows the number of distinct FST states which are stabilized in the equilibrated system. For instance, a "0" signifies that we obtain no FST states (i.e., we have the FM state), while a "4" denotes that all of the 1MML FST, 2MML FST, 3MML FST, and 4MML FST states can be stabilized and distinctly identified. Note that we only record the total number of distinct states stabilized in each case; we do not distinguish between different combinations of FSTs.
According to the results shown in Figure 1b, there is a noticeable transition from the no FST state at a low out-ofplane magnetic field and high DMI constant to four distinct FST states at a high field and low DMI constant. There is a broad parameter window for stabilizing four distinct FST states in the 4MML systems without interlayer exchange coupling. However, with increasing the interlayer exchange coupling, the parameter window for four FST states shrinks, while the parameter windows for three FST states and one/two FST states expand toward the initial 4 FST phase and 0 FST phase, respectively (see also the phase diagram at higher interlayer exchange coupling constant J Interlayer = 0.08 to 0.12 mJ m −2 in Figure S1 of the Supporting Information). We note that even small values of the interlayer exchange coupling are sufficient to substantially alter the stability of FSTs. We can understand these changes by considering the roles of DMI and external magnetic field, which tend to favor perpendicular and parallel spin configurations, respectively. At a given value of the external field, if the DMI is too large, then the labyrinth domain state becomes favorable, and if it is too small then colinear spins are favored. It follows that the effect of the interlayer exchange, which favors parallel spins, is to reduce the required external magnetic field and increase the allowable DMI for FSTs to be stabilized. Hence, the interlayer exchange coupling has the effect of widening the parameter windows for stabilizing 1MML, 2MML, and 3MML FSTs. We additionally observe that for sufficiently large external field and sufficiently small DMI (i.e., the bottom right corner of the phase diagram), the stability of FSTs is diminished by the presence of interlayer exchange coupling. Here, the combined effect of the interlayer exchange and large external field means that larger values of the DMI are required in order to stabilize chiral spin textures such as FSTs. Therefore, by modifying the magnetic parameters of the material system, we can effectively tune the distinct state phase diagram to achieve the maximum number of distinct FST states.
Apart from the results exhibited in Figure 1b, we also performed simulations on an 8MML film to demonstrate the extension of the FSTs in even higher-level cascaded MMLs. As shown in Figure S2 of Supporting Information, we can obtain eight distinct FST states with a large parameter window in the 8MML film. A similar transition to that seen in Figure 1b is also observed. Similarly, by tuning the interlayer FM exchange coupling J Interlayer from 0 to 0.08 mJ m −2 , we observe an apparent shrinking of the region corresponding to the larger number of distinct FST states. In summary, the data in Figure  1b show that to maximize the number of distinct FST states, a relatively high external magnetic field (70 to 100 mT), a moderate DMI constant (1.8 to 2.2 mJ m −2 ), and a relatively weak interlayer FM exchange coupling (lower than 0.08 mJ m −2 ) are required.

Thermal Stability of FSTs.
To reveal the reliability and feasibility of using multiple FSTs in a single MML device, it is important to assess their thermal stability. There exist multiple numerical approaches for characterizing the thermal stability of magnetic textures. The most widely used is the geodesic nudged elastic band method, 38 which can identify the internal energy difference between an initial and a final spin texture. In the present work, we provide a simple analysis of the thermal stability based on the Arrhenius−Neél law by using the difference in internal energy between the intact and decoupled FST states to be their binding energy. By calculating the binding energy, it is possible to characterize the thermal stability of FST states at finite temperatures. We thus determined the binding energy of different FSTs as a function of the in-plane separation between two nearest-neighbor skyrmions in the two topmost layers of the FST stack, respectively. The binding energy of FSTs was calculated according to ref 39; the results are shown in Figure 1c. First, we relaxed the skyrmions in a 2MML FST, 3MML FST, 4MML FST, and 6MML FST, respectively. Then, we artificially shifted the top-layer skyrmion from −50 to 50 nm in the x-direction while fixing the position of skyrmions in the other layers. The total micromagnetic energy of the whole system at every x-shifted position is calculated without relaxing and equilibrating the system. As shown in Figure 1c, we normalized the results by subtracting the energy of initial FST states from the energy after shifting the top layer skyrmion. Therefore, the binding energy as a function of x-shifted position in Figure 1c is greater or equal to 0. The relative value represents the increased total energy induced from shifting the top layer skyrmion.
The binding energy in Figure 1c shows that by shifting the top layer skyrmion away from the center, the total energy first rises gradually and then remains approximately constant after an x-displacement of ∼35 nm. The difference between the highest energy state and the initial state is the binding energy of the FST, which quantifies the energy barrier separating the FST and decoupled skyrmion states. The interlayer exchange coupling J Interlayer = 0.02 mJ m −2 , external magnetic field of 90 mT, and DMI constant of 1.9 mJ m −2 were used when simulating FSTs, as shown in Figure 1c. The binding energy E b of the 2MML FST, 3MML FST, 4MML FST, and 6MML FST is calculated as 1.26, 2.46, 4.60, and 12.81 eV, respectively. The increased binding energy from 2MML FST to 6MML FST can be attributed to enhanced dipolar coupling.
To demonstrate the tunability of the FST thermal stability, we next modified the interlayer exchange coupling constant J Interlayer from 0 to 0.08 mJ m −2 for each FST. The binding energy of FSTs is calculated and displayed in Figure S3 of the Supporting Information. It is observed that the binding energy of each FST follows a linear relationship with the interlayer exchange coupling of the MMLs. To explore this, we calculated the binding energy of a 4MML FST with J Interlayer , ranging from 0 to 0.2 mJ m −2 . The results, as shown in Figure 1d, verify that the binding energy increases with interlayer exchange coupling. Similar results can also be observed for 2MML, 3MML, and 6MML FSTs, as shown in Figure S4 of the Supporting Information.
The larger calculated binding energy results in an improved thermal stability. The thermal stability can be quantified using the Arrhenius−Neél law to estimate the lifetime of a metastable state 38,40 (2) where f = τ 0 −1 is the attempt frequency, k B is the Boltzmann constant, T is the temperature under consideration, and E b is the binding energy, as given by the analysis above. Here, we assume T = 300 K in order to estimate the FST lifetime at room temperature. An attempt frequency of 10 9 −10 12 Hz is typically used. 38,40,41 However, there is some debate about the correct value of the attempt frequency to use, 42 which can be as large as 10 21 Hz. Precisely estimating the lifetime of quasiparticles is beyond the scope of this paper. We therefore choose a large value of 10 21 Hz for the attempt frequency here and give a conservative estimation of the lifetime of FSTs. The binding energy of FSTs extracted from Figure 1c with J Interlayer = 0.02 mJ m −2 leads to estimated lifetimes of 20.1 ns, 1.47 s, 2.12 × 10 20 s, 1.89 × 10 56 s, and 1.58 × 10 194 s for the 1MML FST, 2MML FST, 3MML FST, 4MML FST, and 6MML FST, respectively. Note that the estimated results merely reflect the thermal stability and annihilation probability, rather than the precise lifetime of FSTs in realistic devices. As we are using multiple FSTs in a single device, the least stable one defines the thermal stability of the whole device. However, we can enhance the lifetime of each FST by 5 orders of magnitude by tuning the magnetic parameters and interlayer exchange coupling. Our results here demonstrate the tunability of the thermal stability of FSTs, a key consideration for experimental and commercial device design.

Magnetic and Topological Properties of FSTs.
To utilize FSTs in realistic devices and applications, it is necessary to evaluate their magnetic and topological properties, especially the propagation behavior under applied electric currents. Figure 1 demonstrates that the various MML FSTs have different dimensions in the x−y plane. We thus stabilized a series of FSTs in an 8MML film to demonstrate the MMLdependent diameter of skyrmions.
In the simulations, we artificially equilibrated skyrmions in the select MMLs in order to obtain a 1MML FST, 2MML FST, 3MML FST, 4MML FST, 5MML FST, 6MML FST, 7MML FST, and 8MML FST. The layers are marked as L1 to L8 for the 8MML film in Figure 2a, where the bottom layer is L1, and the top layer is L8, such that the stabilized 1MML FST contains only one skyrmion in L1 while L2 to L8 remain in the saturated FM state. Similarly, the 4MML FST has skyrmions within L1 to L4, and 8MML FST has skyrmions within all layers as a full tube. The relative size and skyrmion diameter of these FSTs can be visualized by superposing 1MML to 8MML FSTs together, as shown in Figure 2a. The skyrmion diameter of FSTs grows monotonically with the number of MMLs. This trend can be attributed to the increased dipolar field as the number of cascaded skyrmions increases. 16 The measured skyrmion diameter of FSTs is shown by the black lines in Figure 2c,d, where an approximately linear relation with the number of MMLs can be observed. The fact that the diameter increases with the number of layers occupied by the FST is likely a result of the enhanced dipolar coupling as the vertical extend of the FST increases.
We then investigated the transport behavior of different MML FSTs under electric current. The spin−orbit torque (SOT) induced by a current perpendicular to the plane (CPP) geometry and the spin-transfer torque (STT) from a current in plane (CIP) geometry are considered. For the CIP case, the charge current is applied in the FM layer in the x-direction with a density of 15 MA cm −2 within each layer of the MML. In the case of CPP, a charge current with density 50 MA cm −2 is applied in the bottom HM layer only; the thickness of the HM layers in each MML is smaller than the electron diffusion length 43 so that no spin Hall effect will be induced. 44 A spin current with spin polarization in the +y direction is created and injected into the first MML starting from the bottom, namely, L1 mentioned above. As a result, FSTs move along a trajectory at an angle to the direction of the applied current, which is the well-known skyrmion Hall effect (SkHE). 45 The angle between the FST trajectory and the direction of applied current (+x in this example) is defined as the skyrmion Hall angle θ SkHE . Under the SOT, FSTs propagate along the +y transverse direction, while under the STT, FSTs move toward the opposite transverse direction. As shown in Figure 2b, the 1MML FST has the largest θ SkHE , while the 4MML FST exhibits a smallest θ SkHE for both CPP and CIP.
To better understand the simulated spin dynamics, it is instructive to employ an analytical model for the motion of non-collinear spin textures. Therefore, we utilize the Thiele equation by imposing the stationary limit that the FSTs move with a constant velocity and that the texture does not deform. We consider both CIP and CPP geometries using the Zhang− Li 46 and Slonczewski 47 torques in the Landau−Lifshitz− Gilbert (LLG) equation (see the Methods section), respectively.
In the CPP geometry, assuming a periodical boundary condition in the x−y plane, the translational motion of spin textures driven by the spin Hall effect can be described by a modified Thiele equation 48,49 ( where G = (0, 0, −4πN) is the Gyroscopic vector with the topological charge N defined in eq 1 and v =(v x ,v y ) is the skyrmion drift velocity along the x and y axes, respectively. The first term G × v in eq 3 is the topological Magnus force that results in the transverse motion of skyrmions as a function of the driving current, which directly results in the SkHE. 45   , and m p is the polarization direction of the spin current. By solving eq 3, we can obtain the skyrmion Hall angle of the FSTs as (4) In CIP geometry, the Thiele equation has the form 11,49 (5) where G, , α, and v are as defined above. β is the nonadiabaticity of STT, for which a small value of 0.02 is used, and j is the current vector describing the direction and amplitude of the applied current. From eq 5, the skyrmion Hall angle of FSTs under CIP is therefore (6) We calculated θ SkHE for both CPP and CIP by eqs 4 and 6, respectively. Due to the fact that different FST states have different diameters, the calculated value of D xx is different for each FST state. Hence, the skyrmion Hall angle varies with the number of magnetic layers occupied by the FST. Comparison of the theoretical prediction and simulation results of θ SkHE for 1MML to 8MML FSTs is exhibited in Figure 2c,d, where Figure 2c shows CPP geometry and Figure 2d shows CIP. The simulated results of θ SkHE are shown by the discrete data points, and the calculated results from Thiele's equation are exhibited as solid lines. There is a decrease of θ SkHE both in simulations and theoretical calculations when cascading MMLs for FSTs in both CPP and CIP geometry. However, there is a slight disparity between the theoretical predictions and simulations for θ SkHE , which becomes more significant as the number of MML repeats increases. The overestimation of skyrmion Hall angle by the Thiele equation can be explained by additional dissipation mechanisms related to dynamic variation of the skyrmion shape, which cannot be captured by the rigid shape approximation used to derive the Thiele equation.
We also extracted the velocity of FST propagation in CPP and CIP from simulations. The results are shown in Figure  2c,d. FSTs generally move faster in CPP than that in CIP, even with a smaller amplitude of the applied electric current. This trend is consistent with previous studies comparing CPP and CIP motion and likely results from the more significant Slonczewski in-plane torque in CPP geometry compared to the smaller, field-like torque which results from CIP. 2,13 On the other hand, a marked drop of FST velocity from 7.5 m s −1 for 1MML FST to 3 m s −1 for 6MML FST under the same driving current can be observed in the CPP geometry, while the velocity of FSTs remains virtually constant as the number of MML repeats increases for the case of CIP geometry. The different current-dependent characteristics of skyrmion velocity could result from the way that torque is injected into the magnetic textures via the CPP and CIP geometry. In the CIP geometry, the torques act on individual skyrmions in each layer, while in the CPP geometry, the torque is merely applied to the top layer skyrmion so that velocity would decrease as the number of MML repeats increases. Overall, these results indicate that the proposed FSTs have tuneable thermal stability and distinct topological and magnetic properties, which may be used as signatures to identify them from each other. In the next chapter, we demonstrate one of the potential uses for FSTs: an MML multiplexing device.

MML Multiplexing Device with FSTs. 2.3.1. Demonstration of the Proposed Multiplexing Device.
In this section, a multilayer multiplexing device is proposed using multiple FSTs as information carriers. Such a device can perform signal multiplexing, signal transmission, and automatic signal demultiplexing. We have shown that as many as eight distinct FST states, possibly more, can be achieved through careful selection of the system and magnetic parameters. In this section, an FST-based device built from a 4MML nanotrack is illustrated for simplicity. The schematic of the proposed 4MML multiplexing device is illustrated in Figure 3. The multiplexing device has three parts: a 4MML nanotrack for FST transmission, terrace-like MML stages for FST nucleation via magnetic tunnel junctions (MTJs) fabricated on the top, and a four-branched track for FST automatic selecting.
Experimental fabrication of this device will require a complex, multiple-step patterning process. For example, the terrace structure for FST nucleation may be fabricated via several subtractive lithography steps, with single-layer removal achieved through carefully calibrated ion-beam milling or focused ion beam etching. Given the modular nature of the device operation, it is in principle possible to optimize each component of its functionality (i.e., the FST stability, nucleation, chambering, propagation, and detection) individually, with the full device operation verified thereafter. Indeed, some of the required steps have recently begun to be experimentally explored, including current-induced skyrmion nucleation in MTJs 50 and voltage-controlled-magnetic-anisotropy (VCMA) control of their transport through nanotracks. 51 The proposed device shown in Figure 3 can multiplex and transmit four distinct sequences of information signals simultaneously, where each information signal is encoded by one type of FST, and the presence/absence of the FST encodes information "1"/"0". As shown in Figure 3, each type Single-layer precision may be achieved either through careful calibration of the material growth/removal rate or by dose-modulated electron-beam lithography. 52 In order to contact the top electrode of the MTJs for tunnelling in the perpendicular direction, each of them must be sheathed by an insulating material such as SiO 2 .
The simulated system has a four-level terrace geometry with cell size 2 nm × 2 nm × 1 nm. From left to right, the geometry of each terrace is 600 nm × 600 nm × 2 nm, 600 nm × 600 nm × 4 nm, 600 nm × 600 nm × 6 nm, and 600 nm × 600 nm × 8 nm. The MML has a small interlayer exchange coupling of J Interlayer = 0.02 mJ m −2 . The simulations were performed at room temperature by including an extra thermal fluctuation field 31 and with room-temperature magnetic parameters. The magnetic parameters were retrieved from the literature of experimental measurements (see the Methods section). Electric current is injected into the MMLs from the MTJ using a CPP geometry, and the FSTs are nucleated through the STT. The pulse width of the injected current is 1.5 ns FWHM with a peak current density of 250 MA cm −2 . Although the pulse length is very short, the requisite current density is very large and may cause damage to the ultrathin tunnel barrier of the MTJ through Joule heating. 53 Therefore, future work is required to minimize the effects of Joule heating, either through optimization of the spin injection process or by including a heatsink.   Figure 4d verifies the results in Figure 4b.
Since a stochastic thermal field was included in the simulations of FST nucleation, the successful nucleation of FSTs occurs probabilistically, rather than being deterministic. We therefore determined the probability phase diagram of the FSTs by scanning the electric pulse width from 0.1 to 12.8 ns in a geometric sequence (i.e., 0.1, 0.2, 0.4, 0.8, 1.6, 3.2, 6.4, and 12.8 ns) and the current density from 100 to 800 MA cm −2 in 50 MA cm −2 increments. We simulated the nucleation process for each FST with each given pulse width and current density 20 times, then calculated the probability of successful nucleation. The probability phase diagrams of successful nucleation of 1MML to 4MML FSTs are shown in Figure  S5 in Supporting Information. There are 15 × 8 = 120 data points displayed in the figures, where each data point illustrates the probability of successful nucleation out of 20 distinct attempts. The color coding describes the probability value, where white represents 0% of the probability and black represents 100%. A transition from low nucleation probability for 1MML FST to high nucleation probability for 4MML FST can be seen in Figure S5. However, for the multiplexing device proposed in this work, four different FSTs need to be nucleated at the four-stage terraces individually at the same time. Therefore, we also calculated the probability for successfully nucleating four FSTs simultaneously, with results summarized in the phase diagram displayed in Figure 4c. The results suggest that FSTs are more likely to be nucleated when we apply a large amplitude of current density (200 to 300 MA cm −2 ) with a relatively large pulse duration (greater than 1 ns). This information is critical for experimental realization of the multiplexer device, wherein controlled generation of individual FSTs is required. There has been much recent progress in experimentally realizing single-skyrmion nucleation in multilayer systems, including via defects generated by ionirradiation, 54 through local applied magnetic fields from a magnetic tip, 55 and through spin-polarized currents. 56 We neglected the spin memory loss (SML) 46 when injecting STT to magnetizations in cascaded MMLs to simplify the nucleation procedure. Careful consideration of SML in the device-level simulations may be the focus of follow-up work. The results in this section demonstrate the nucleation of FSTs in the proposed multiplexing device, including the structure of the four-stage terrace nucleation site, the detailed nucleation process of FSTs, and the probability phase diagram of successful nucleation. In the following section, the second important stage of the device workflow will be addressed, namely, the FST initialization process.

Chambering the Multiplexing Device with FSTs.
After nucleating FSTs in the four-stage terraces, the proposed device needs to be initialized before all of the FSTs propagate in the main 4MML nanotrack. The FSTs must propagate upward from the nucleation sites into the main nanotrack during the initialization process. We refer to the initialization process as "chambering" in this work. FSTs must overcome an energy barrier when moving into the 4MML nanotrack from a film with a different number of MMLs, 39 which can be attributed to the dipolar field in the two systems. As shown in Figure 5a, the 4MML multiplexing device needs to be chambered with all four FSTs from the four-stage terraces. We used CIP geometry for the chambering process because STT can act on every skyrmion throughout the FSTs, which decreases the risk of FST decoupling. Figure 5a illustrates the chambering process enabled by an STT resulting from the CIP in the y-direction with the current density of 100 MA cm −2 and pulse width of 4 ns. FSTs first propagate toward the boundary of different MML systems in the y-direction, i.e., along the direction of the applied current, before moving along the boundary in the x-direction for a short distance. They subsequently cross the boundary and successfully move into the main track. The 1MML FST exhibits the largest displacement along the x-direction; this displacement gradually reduces as the number of MML repeats increases. We then varied the interlayer exchange coupling constant J Interlayer and the amplitude of the current density to generate a phase diagram of the initialization process. By changing J Interlayer and current density, the initial state leads to four phases which we label "chamber", "stuck", "nucleate", and "decouple", as shown in the left panel of Figure 5b. If the interface traps the FST, we mark it as a "stuck"; if all four FSTs propagate into the main track, we mark it as a "chamber"; if any of FSTs annihilate, we mark it as a "decouple"; if the skyrmions in different layers separate, we also mark this as a "decouple"; and if a new skyrmion is nucleated within the FM layers during the process, we mark it as a "nucleate".
The results of the initialization process after injecting electric current are summarized in the phase diagram shown in Figure  5b. The four possible phases "chamber", "stuck", "nucleate", and "decouple" are marked with a blue rectangle, a red triangular, a green circle, and a yellow star, respectively. As the target phase, the "chamber" phase has a decently large parameter window, which occupies the upper left 1/4 of the phase diagram. In the multilayers with large interlayer exchange coupling constants (J Interlayer ≥ 0.06 mJ m −2 ), 1MML FSTs are stuck by the boundary when applying insufficient current densities (J dc < 50 MA cm −2 ) and will annihilate when the current density is large (J dc ≥ 50 MA cm −2 ). 3MML FSTs will nucleate a new skyrmion in the top FM layer when crossing the boundary with high interlayer exchange coupling constants (J Interlayer > 0.03 mJ m −2 ) and applied current densities (J dc > 37.5 MA cm −2 ). 2MML FSTs behave more robustly and produce the "chamber" phase in most situations except for at extremely high interlayer exchange coupling constants (J Interlayer > 0.08 mJ m −2 ) and current densities (J dc > 125 MA cm −2 ), wherein annihilation is observed. 4MML FSTs always produce the "chamber" phase, resulting from the fact that there is no energy barrier during the propagation process in this case. The superposed markings in the phase diagram indicate different behaviors from different FST types. For instance, all cases of the yellow stars superposing green circles correspond to the situation that the 1MML FST annihilates, and the 3MML FST nucleates a new skyrmion; the red triangle with the green circle indicates that the 1MML FST is stuck, and the 3MML FST nucleates a new skyrmion.
In order to better explain the phases shown in Figure 5b, we extracted the total micromagnetic energy of the regions containing 1MML FST, 2MML FST, 3MML FST, and 4MML FST individually. The trajectories when chambering 1MML FST, 2MML FST, and 3MML FST are marked as C1, C2, and C3, respectively, as shown in Figure 5a. The energy evolution of C1, C2, and C3 can be seen from Figures 5c, S6 of Supporting Information and Figure 5d. Here, we varied the interlayer exchange coupling constant from 0 to 0.08 mJ m −2 while fixing the current density J dc = 125 MA cm −2 . With a higher interlayer exchange coupling constant, the "decouple" phase is observed in C1, while the "nucleate" phase could be seen in C3. The energy evolution is further demonstrated by the change of topological charge of the system, as shown in the insets of Figure 5c,d. The energy barrier for chambering 1MML FST into the central track increases as a linear relationship with the amplitude of J Interlayer . When chambering 1MML FST, the "chamber" phase is obtained when J Interlayer ≤ 0.06 mJ m −2 ; therefore, no net change of topological charge is observed. Annihilation of the skyrmion happens at J Interlayer = 0.08 mJ m −2 , so the topological charge of the system changes to zero. Similarly, when chambering the 3MML FST, a net increase of the topological charge by around 1/3 occurs, suggesting the nucleation of a new skyrmion in the system at higher interlayer exchange coupling constants (J Interlayer ≥ 0.04 mJ m −2 ). Note that the initialization process of 2MML FST, marked as C2 in Figure 5a, merely exhibits "chamber" and "stuck" phases in Figure 5b, and therefore, no noticeable variation can be observed in the topological charge. The energy evolution when chambering the 2MML FST is presented in Figure S6 of the Supporting Information.

Voltage-Controlled Synchronizers for Pipelining.
After the initialization process, we now have the multiplexing device chambered with FSTs in the main nanotrack. FSTs can then be transmitted along the track toward the demultiplexing region. Considering the long transmission distances likely required in real-world applications, it is better to divide the transmission track into several regions, as shown in Figure 6a. In such a design, we can achieve pipelined transmission for information carriers to enhance the device throughput. The results in Figure 2 have demonstrated that FSTs with different MMLs propagate with various velocities, where the 4MML FST moves the most slowly, and the 1MML FST the most quickly. Therefore, synchronizers are required to maintain the correct order of information sequences and avoid data loss during transmission. 25 Here, we can utilize voltage gates as synchronizers based on the VCMA effect. The VCMA effect was first reported in a 3D transition FM layer in 2−4 nm thick FePt (Pd) films. 57 Surprisingly, it has been reported that a small electric field of 100 mV nm −1 is sufficient to change the PMA by 40%, corresponding to a VCMA efficiency of 210 fJ V −1 m −1 at room temperature. 58 In this work, the simulation of the VCMA effect is based on a linear relationship 12 (7) where ϑ is the VCMA coefficient, V b is the bias voltage on the VCMA gate, and K u0 is the background anisotropy constant. Regions in which the PMA is elevated provide an increased energy barrier, whereas regions of reduced PMA lead to a potential well. When simulating the pipelined transmission of FSTs schematized in Figure 6a, we periodically set VCMA gates on and off by changing the amplitude of the PMA constant by 10%. An electric driving current with the amplitude of 3 MA cm −2 is applied to the HM layer in the x-direction. Compared to the chambering process of FSTs where STT is utilized, we apply SOT for the transmission of FSTs instead to achieve higher driving efficiency and thus lower energy consumption. 13 The pipelined transmission of four FSTs in the 4MML main track is illustrated in Figure 6a. We deployed three VCMA gates in the track, which are marked as red shaded rectangles. The track is divided into four regions as registers, each of which can store a packet of information consisting of four bits. The presence (absence) of a 1MML FST/2MML FST/3MML FST/4MML FST encodes a single binary "1" ("0") for the first/second/third/fourth bit, respectively. Initially, four packets of information signals "1111", "1110", "0011", and "0001" are positioned in consecutive registers of the device via propagation mediated by the VCMA gates. They are then transmitted along the x-direction and synchronized by adjacent VCMA gates. It should be noted that the FSTs are nucleated and propagated from left to right in order of the number of MML repeats, such that within an information packet, the order is [1MML FST|2MML FST|3MML FST|4 MML FST]. This design can retain the order of FSTs and prevent mismatch between information packets during transmission. Figure 6a indicates that the information integrity and sequence order are well retained, allowing the transmission efficiency to be quadrupled. Note that for clarity, Figure 6a shows a schematic representation of the VCMA-controlled FST transmission. The full simulation results are provided in Video S1 of the Supporting Information. Furthermore, this pipelined design is not limited to four-bit packets of information. The VCMA gates can support as many information carriers in the nanotrack as can be stabilized in the MML system. Here, we use four carriers in each packet for illustration and consistency with the other results in this paper.

Automatic Demultiplexing via Track Selectors.
The final stage of the proposed MML multiplexing device workflow is to filter the FSTs out of the information packets to decode the information signals. Here, we propose to use a four-branch track selector to filter the FSTs into individual tracks for detection. We can utilize a track selector to filter different FSTs, as illustrated in Figure 6b. There are two choices to filter the FSTs: the first one is via the skyrmion Hall angle θ SkHE and the second is the size difference of the skyrmion diameter. Changing the angle or the width of the branch can filter and demultiplex different FSTs. As for the skyrmion Hall angle θ SkHE , our results suggest that there is a difference in θ SkHE for FSTs with different MMLs. Namely, θ SkHE for 1MML FST to 4MML FST under SOT is −47.5, −44.0, −40.7, and −35.0°, respectively. At the same time, there is also a difference in the skyrmion diameter among FSTs, where the skyrmion diameter for 1MML FST to 4MML FST is 21.5, 26.5, 32.5, and 38.0 nm, respectively. The angle-based selector would require a much larger device footprint to filter the FSTs for their relatively small angle difference. Therefore, in this work, we propose to demultiplex the information signals by a widthbased track selector.
Simulations of the demultiplexing process were performed to obtain a phase diagram for the four-branch track selector. Here, we modified the width of the branch from 50 to 130 nm while the angle of the branch was fixed at 70°, and the electric current density was 3 MA cm −2 (see Figure 6b). The four FSTs are individually propagated through the selector along the track, driven by the SOT. If the FST propagated into the branch, we marked the result as "go upward". Otherwise, it was marked as "go straight". The track selecting phase diagram for 1MML to 4MML FSTs is shown in Figure 6d, where the skyblue colored region represents FSTs propagating into the branch and the cream-yellow colored region denotes FSTs moving straight, passing by the branch. The results indicate enough space to design a track selector for filtering four FSTs out of one packet because each FST has an exclusive window with 20 nm width. Guided by the track selecting phase diagram for FSTs, we configured a four-branch track selector, as shown in Figure 6c. The width of these four branches is 60, 80, 100, and 120 nm for demultiplexing 1MML FSTs, 2MML FSTs, 3MML FSTs, and 4MML FSTs, respectively. Successful demultiplexing of the FSTs, as shown in Figure 6c, verifies the feasibility and validity of such a width-based track selector.
In summary, these results demonstrate that the proposed FSTs in the MML system can be potential candidates as information carriers, based on the fact that multiple FSTs can be nucleated, transmitted, and filtered in a single MML device. The results highlight the potential for three-dimensional skyrmionics quasiparticles to significantly expand the information storage and processing capability of tailored MML systems.

CONCLUSIONS
This investigation was designed to assess the hypothesis that in MML systems, skyrmions can exist within part of a multilayer as FSTs. We confirmed this with magnetic energy analysis and micromagnetic simulations. The findings suggest that distinct FST states may coexist in a single MML system which are tuneable in their thermal stability and magnetic properties. Their topological properties and current-driven behavior were analyzed both with theoretical calculations and micromagnetic simulations. This work also proposes to use such FSTs to encode information in an MML multiplexing device, where multiple FSTs can be nucleated, transmitted, and filtered. This study provides the first comprehensive assessment of encoding information by multiple FSTs in a single MML device, highlighting the potential utility of distinct skyrmion states in MML systems. Further investigations are required to explore device settings for FSTs and establish effective nucleation and detection methods.

Micromagnetic Simulations.
The micromagnetic simulations were performed using the GPU-accelerated micromagnetic program Mumax 3 . 31 The time-dependent magnetization dynamics are conducted by the LLG equation (8) where m = M/M s is the reduced magnetization, M s is the saturation magnetization, γ LL is the gyromagnetic ratio, h eff = H eff /M s is the reduced effective field, α is the damping parameter, is the SOT efficiency with T −1 s −1 being the gyromagnetic ratio of an electron, ℏ is the reduced Planck constant, j e is the current density, θ SH is the spin Hall angle, e is the electron charge, t is the thickness of the FM layer, and m p is the polarization direction of the spin current. The energy density E is a function of m, which contains the exchange energy term, the anisotropy energy term, the Zeeman energy term, the magnetostatic energy term, and the DMI energy term. The material parameters to perform the simulations are chosen according to previous reported roomtemperature experimental results: 20 damping parameter α = 0.1, DMI constant D int = 1.5 mJ m −2 to 2.6 mJ m −2 , the value for the Gilbert gyromagnetic ratio γ = −2.211 × 10 5 mA −1 s −1 , saturation magnetization M s = 956 kA m −1 , the spin Hall polarization Θ SH = 0.6 to enhance the spin Hall effect, the uniaxial out-of-plane magnetic anisotropy K u = 717 kJ m −3 , the polarization of the spin current is in the +y direction, and the exchange constant is assumed to be A = 10 pJ m −1 . To ensure the accuracy of calculation, the mesh size is set to 2 nm × 2 nm × 1 nm, which is smaller than the exchange length nm and DMI length nm. An external magnetic field of 10 to 100 mT in the out-of-plane direction is applied for the simulations. In the simulated MMLs, the intermediate HM 1 and HM 2 layers are thinner than the electron spin diffusion length. In this case, the torques would be efficient only in the external layers. In the simulation of FST propagation along the nanotrack, the SOT created via a CPP is applied only in the bottom layer, and the injected spin polarization is uniform in the layer. The injected current is then modeled as a fully polarized vertical spin current. In the simulation of the chambering process, electric current is applied in the +y direction via the CIP geometry, where skyrmions in all FM layers within FSTs are driven by STT.