Forming individual magnetic biskyrmions by merging two skyrmions in a centrosymmetric nanodisk
Börge Göbel Jürgen Henk Ingrid Mertig


[
    [
    December 30, 2023
=====================

empty



§ INTRODUCTION


Magnetic skyrmions <cit.> are whirl-like magnetic quasiparticles on the sub-micrometer scale, first observed in MnSi <cit.> as periodic arrays (so-called skyrmion crystals or lattices). 
Later, also isolated skyrmions have been found in ferromagnetic Fe_0.5Co_0.5Si films <cit.>. Each skyrmion is characterized by a topological charge 

    N_Sk= 1/4π∫m⃗(r⃗) ·[ ∂m⃗(r⃗)/∂ x×∂m⃗(r⃗)/∂ y] dr^2

of ± 1 [m⃗(r⃗) is the magnetization], which imposes an energy barrier that protects a skyrmion from annihilating to the ferromagnetic groundstate. Their high stability and small sizes make skyrmions candidates for the carriers of information in future storage devices <cit.>. For example, skyrmions can be written and deleted, driven by electric currents and read in thin films <cit.>.

An object which is closely related to the skyrmion is the biskyrmion, a term that describes two skyrmions in a bound state. Like the skyrmion <cit.> it was initially proposed in nuclear physics <cit.>. In the context of magnetism the term `biskyrmion' has been used for two skyrmions in a bilayer quantum Hall system <cit.>, for two asymmetric skyrmions in the cone phase of a chiral helimagnet which exhibit an attractive interaction <cit.>, and – as observed for the first time in 2014 <cit.> – a composition of two partially overlapping skyrmions with opposite in-plane magnetizations (Fig. <ref>b) <cit.>. Still, all of these objects are geometrically distinct. In the remainder of this paper the term `biskyrmion' refers always to the latter object: a pair of partially overlapping skyrmions with a helicity difference of π that has attracted an enormous research interest since its initial discovery. 

Although merged, each skyrmion forming the biskyrmion can still be identified and their shifted centers indicate a merely partial overlap. For this reason, each skyrmion contributes with its topological charge of ± 1 to the biskyrmion's topological charge of N_Sk = ± 2. Periodic arrays of such biskyrmions (biskyrmion crystals) have first been observed in the centrosymmetric layered manganite La_2-2xSr_1+2xMn_2O_7 <cit.>; other centrosymmetric biskyrmion hosts followed recently <cit.>. 



If individual biskyrmions could be stabilized, they could be used to store information just like conventional skyrmions. However, up to now, individual biskyrmions have only been predicted in frustrated magnets <cit.> — even the existence of conventional skyrmions in these materials <cit.> remains to be proven experimentally. These considerations suggest to investigate centrosymmetric materials with dominating dipole-dipole interaction, like in the experimental observations of biskyrmion crystals, to find also individual biskyrmions. In these materials the constituents of biskyrmions — skyrmionic bubbles — have been found already decades ago <cit.>.

In this Paper, we propose a scheme for creating biskyrmions in a controlled manner. In a first step, we derive an attractive skyrmion-skyrmion interaction which is necessary for the formation of a biskyrmion, using an analytical model. Thereafter, we show by micromagnetic simulations for a nanodisk geometry that two Bloch skyrmions with opposite in-plane magnetizations can be written at opposite sides of a nanodisk. We simulate the resultant motion of such skyrmions upon relaxation and derive their equations of motion analytically. After relaxation we find that both skyrmions have merged to a stable biskyrmion state (Fig. <ref>). The importance of the skyrmions' helicities for this process is discussed. We summarize our findings and give an outlook.




§ RESULTS




 §.§ Analytical superposition of two skyrmions

As presented in the introduction, it is an established fact that non-collinear spin textures like skyrmions can be stabilized by the dipole-dipole interaction <cit.>. Therefore, it is well conceivable that also two skyrmions with reversed in-plane magnetizations can be stabilized in the same sample by this mechanism. The question is now how these skyrmions interact. As we will show, already considering a short-range approximation of the dipole-dipole interaction leads to the conclusion that the two skyrmions attract each other, thereby allowing to form a biskyrmion. This finding is later confirmed by micromagnetic simulations where the full dipole-dipole interaction is considered and a biskyrmion state is metastabilized. 

To establish the attractive interaction between two skyrmions with opposite in-plane magnetizations, we superpose such two (fixed) skyrmions at a distance of Δ r analytically. For simplicity we neglect shape deformations of the skyrmions (such deformations are accounted for by the micromagnetic simulations later in this paper). 
The magnetic textures  m⃗(r⃗) = (sinΘcosΦ, sinΘsinΦ, cosΘ) are imposed onto a two-dimensional square lattice with lattice constant a. 

The spherical coordinates of the magnetic moments of an individual skyrmion read Φ = arctan(m   y/x) + γ and Θ = 2 arctan(√(x^2 + y^2) / r_0); r_0 determines the size of the skyrmion. The helicity γ of the skyrmion manifests itself in the in-plane components via Φ = m ϕ + γ. Here, ϕ is the polar angle of the position vector with respect to the skyrmion's center and m = +1 is the vorticity, which defines the sense of in-plane spin rotation. In case of a Bloch skyrmion we have m=1 and γ = ±π / 2, i. e., the polar angle of the spin orientation Φ has a fixed offset to the polar angle ϕ of the position r⃗ with respect to the skyrmion's center. The topological charge N_Sk = p · m = ± 1 of a skyrmion depends on the skyrmion's polarity p=± 1. In this Paper, we consider skyrmions in a ferromagnetic surrounding, magnetized along -z, which gives a polarity and topological charge of  +1.

One way to superpose two skyrmions is adding their polar angles and multiplying the arguments of their arctan function for the azimuthal angle, 

    Φ    = arctan(y-Δ y/2/x-Δ x/2)+arctan(y+Δ y/2/x+Δ x/2)+(γ_d-α), 
    Θ    = 2arctan[1/r_0^2√((x-Δ x/2)^2+(y-Δ y /2)^2)·√((x+Δ x/2)^2+(y+Δ y /2)^2)].

The condition |m⃗| = 1 holds for all r⃗. Here, Δr⃗ = (Δ x, Δ y, 0) = Δ r (cosα, sinα, 0) is the displacement vector of the two skyrmions with helicities γ_d and γ_d+π. For an inter-skyrmion distance of Δ r ≫ r_0 we find two isolated skyrmions, for Δ r ≈ 2r_0 the solution is a biskyrmion state with two skyrmions overlapping partially, and for Δ r = 0 a higher-order skyrmion is formed. In any case, the total topological charge is N_Sk = +2.



 §.§ Attractive skyrmion-skyrmion interaction

We clarify now which term in the lattice Hamiltonian causes an attractive interaction of the two skyrmions, what can lead to the formation of a biskyrmion. The exchange interaction H_ex = -1/2∑_ij J_ijm⃗_i ·m⃗_j is independent of Δ r and independent of the skyrmions' sizes and shapes. Hence, the exchange interaction alone will not stabilize a biskyrmion. This finding is in line with the fact that even conventional skyrmions cannot be stabilized by a bare exchange interaction, unless further (frustrated) exchange interactions are considered; then, the interaction constants J_ij have different signs for nearest and second-nearest neighbor spins. The same requirement holds for biskyrmions<cit.>.

In most materials, skyrmions are stabilized by the Dzyaloshinskii-Moriya interaction (DMI) <cit.> H_DMI = 1/2∑_ijD⃗_ij· (m⃗_i ×m⃗_j) that originates from spin-orbit coupling and a broken inversion symmetry. The DMI vectors D⃗_ij = -D⃗_ji whose directions are prescribed by the crystal symmetry determine the metastable spin textures <cit.>. For all types of DMI the energy contributions of the two skyrmions with reversed in-plane magnetizations — what is necessary for the formation of the here considered biskyrmion — cancel. As a result, H_DMI is independent of Δ r and does not allow for the formation of a biskyrmion. Even more detrimental, at least one of the two skyrmions would annihilate in a simulation.



An interaction that is present in all magnetic materials, but `overshadowed' by the DMI in non-centrosymmetric materials, is the dipole-dipole interaction

    H_dp    = -1/2∑_ijμ_0/4π |r⃗_ij|^5[3(m⃗_i·r⃗_ij)(m⃗_j·r⃗_ij)-m⃗_i·m⃗_j].
 
In case the DMI is forbidden, like in centrosymmetric materials, this interaction can stabilize not only skyrmions <cit.> but also biskyrmions, as will be presented later in this paper by micromagnetic simulations. To establish an intuitive understanding of the attractive skyrmion-skyrmion interaction, which is necessary for the formation of a biskyrmion, we consider a short-range approximation of the dipole-dipole interaction for the calculations in this section. Note however, that this approximation is accompanied by a loss of information <cit.>: Even the lowest energy state of the considered model would relax to a collinear ferromagnetic configuration if only the short-range approximation was considered in a simulation in which shape relaxation is allowed. Since Derrick's theorem <cit.> implies the absence of local energy minima for non-collinear spin textures for this reduced type of dipole-dipole interaction <cit.>, we cannot make statements on the metastability of a certain configuration. However, comparing the energies of states described by the model in Eq. (<ref>) allows to extract an attractive skyrmion-skyrmion interaction. The actual proof of metastability are the micromagnetic simulations of the three-dimensional system with the full dipole-dipole interaction shown later in this paper. A three-dimensional consideration is essential since non-collinear spin textures are only metastable once the sample is thicker than a critical value <cit.>. 

To extract the attractive interaction we consider in this section only nearest-neighbor interactions in a two-dimensional square lattice. The second term in Eq. (<ref>) merely rescales the exchange interaction, thereby rendering only the first term relevant

    H_dp,approx∝-∑_⟨ij|⟩(m⃗_i·r̂⃗̂_ij)(m⃗_j·r̂⃗̂_ij),

where r̂⃗̂_ij is a unit vector which is parallel to either x̂⃗̂ or ŷ⃗̂. This interaction can be interpreted as a special easy-plane anisotropy. It favors parallel spins energetically. Such an alignment is approximately fulfilled for the overlap region of the two skyrmions (that is the center of the biskyrmion). As a result, two skyrmions in centrosymmetric materials attract each other.

We find that partially overlapping Bloch skyrmions (γ=±π/2) at an displacement angle of α=±π/4 or ± 3π/4 have the lowest energy of all states described by Eq. (<ref>) with respect to the short-range approximation of the dipole-dipole interaction. The energy minimum occurs at an inter-skyrmion distance of Δ r ≈ 2   r_0, which is a biskyrmion state (Fig. <ref>a). The helicity dependence (Fig. <ref>b) has a cosine shape making the Bloch skyrmion energetically most favorable. This helicity dependence is not unique to the biskyrmion state; a similar dependence is found for a single skyrmion. The α dependence (Fig. <ref>c; global rotation of the texture) originates from the underlying lattice. As a consequence it has a period of π/2 with respect to α for the square lattice. Different lattices have energy minima at other angles, e. g. α = ±π/6, ± 3π/6 and ± 5π/6 for a hexagonal lattice.

So far, we have shown that two skyrmions with reversed in-plane magnetizations attract each other in centrosymmetric systems without DMI but with dipole-dipole interaction. The resulting biskyrmion is expected to be of Bloch type and to form at an angle of 45^∘ with respect to the underlying lattice. In the following, micromagnetic simulations will be used to analyze how these two specific skyrmions can be stabilized in a three-dimensional sample considering the full dipole-dipole interaction. We will show that all three results from the analytic model are confirmed in the simulation: the skyrmions attract each other to form a biskyrmion characterized by γ_e=±π/2 and α=±π/4,± 3π/4.



 §.§ Setup for the micromagnetic simulations

To support the above analytical considerations, we performed micromagnetic simulations for a thin magnetic nanodisk with out-of-plane anisotropy. These disks have been in the focus of skyrmion-related research in many recent publications <cit.>. The considered disk and its surfaces exhibit negligibly small DMI, so that only the dipole-dipole interaction (included in the demagnetization field) could stabilize a biskyrmion. 

The micromagnetic simulations were performed using the Mumax3 code <cit.>. Therein, the Landau-Lifshitz-Gilbert equation (LLG) <cit.> with parameters similar to those in Ref. <cit.> is solved for each discretized magnetic moment m⃗_i; see Methods.

Initially, two energetically degenerate metastable skyrmions with helicities ±π/2 are `written' symmetrically to the nanodisk's center. Upon relaxation they are pushed to the center by their attractive interaction and the confinement potential of the nanodisk. When the skyrmions reach the energetically optimal distance, a biskyrmion has formed. We proceed by discussing the individual steps in detail and begin with the mechanism for writing the skyrmions.



 §.§ Generation of isolated skyrmions

The vanishing DMI complicates a controlled writing of skyrmions in centrosymmetric materials, since there is no energetically preferred chirality: injecting spins oriented along z⃗ via spin torque, as commonly done in DMI-dominated systems <cit.>, produces a topologically trivial bubble (the in-plane spin components point toward the disk's center; they are oriented nearly parallel). Evidently, in order to write two skyrmions off-center without DMI the aimed-for texture has to be induced explicitly. In other words, deterministic instead of stochastic writing has to be considered. This is especially important since the two skyrmions need to have opposite helicities. 

We consider two smaller disks attached to the nanodisk (Fig. <ref>a). These additional disks host stable in-plane vortex configurations <cit.>, similar to Refs. <cit.>. The opposite helicities of the vortices, which are switchable e. g., by a magnetic field <cit.>, determine the helicity of the two emerging Bloch skyrmions. If an electric current j⃗_write is applied in perpendicular direction (here: j⃗_write =7× 10^8 A/cm^-2e⃗_z applied for 200 ps), the spin of the flowing electrons is aligned with the vortex (as long as the coupling is sufficiently large) and injected into the underlying nanodisk; there a skyrmion is formed via spin torque.



The writing process is visualized in Figs. <ref>b–g. Starting from a ferromagnetic state, a skyrmionium-like configuration (a skyrmion with reversed skyrmion-like magnetization in its center <cit.>) is written (b–d). Since the `central' reversed skyrmion is very small, it quickly annihilates (< 100 ps) so that a skyrmion with a topological charge of N_Sk=+1 remains at each side of the disk (e). On the time scale of 1 ns their sizes relax (f, g).

While the writing mechanism considered here appears convenient, other approaches are suited for an experimental realization as well; in the following, we give an incomplete list.

Merons (vortices with out-of-plane spins in their center) in the writing devices work equally well. The same holds for skyrmions. All of the considered textures have a similar in-plane magnetization which seems to be the most important component for skyrmion generation in centrosymmetric materials; due to the fact that the dipole-dipole interaction is achiral, a suitable in-plane magnetization has to be `imprinted' in order to `generate' a topological charge.

In a recent publication <cit.> we have presented a mechanism for writing skyrmioniums via a photosensitive switch <cit.>: a laser triggers a radial current in a ring of a heavy metal, thereby generating an in-plane toroidal spin polarization via the spin Hall effect. This toroidal spin profile is then injected into the magnetic nanodisk, very similar to the vortex-writing device but for a ring geometry instead of a disk geometry. Depending on the chosen parameters, this method works for writing skyrmions as well. The sign of the applied bias voltage controls the spin polarization's orientation and thus the writing of either +π/2 or -π/2 Bloch skyrmions.

As a third option, one may consider spatiotemporally tuned electron sources <cit.>. The type of induced Bloch skyrmion can be controlled via the propagation direction of the electrons (beam applied from above or below the nanodisk).



 §.§ Motion of a single skyrmion

Before discussing the formation of biskyrmions, we address the motion of a single skyrmion within the nanodisk. After relaxation of shape and size, a skyrmion starts to move toward the center of the nanodisk in a counter-clockwise spiral trajectory (orange in Figs. <ref>a–c). Initially, the skyrmion is accelerated but is subsequently decelerated when the skyrmion's center is close to the nanodisk's center. The helicity of the skyrmion remains unchanged during the evolution.



The spiral trajectory is understood by considering the Thiele equation <cit.> 
    b G⃗×v⃗ - b Dαv⃗ + F_int(r) r⃗/r = 0.

In this effective `center-of-mass-like' description of magnetic quasiparticles with velocity v⃗ all internal degrees of freedom have been integrated out. The spatial details of the skyrmion are `condensed' in the gyromagnetic coupling vector G⃗ = -4π N_Ske⃗_z and in the dissipative tensor D with elements D_ij = ∫∂_is⃗(r⃗)·∂_js⃗(r⃗) d^2r. Only  D_xx = D_yy are nonzero. b=M_sd_z/γ_e is determined by the saturation magnetization M_s, the disk thickness d_z and the gyromagnetic ratio of an electron γ_e. Due to the disk geometry the force F_int, that covers interactions of the considered quasiparticle with other quasiparticles or with the edge of the sample, is radially symmetric. The Thiele equation in polar coordinates

    4π N_Skb(rϕ̇e⃗_r-ṙe⃗_ϕ)-α D_xxb(ṙe⃗_r+rϕ̇e⃗_ϕ)+F_int(r)e⃗_r = 0

 yields 

    ṙ   =F_int(r)/(α D_xx)^2b+(4π N_Sk)^2bα D_xx,
    ϕ̇   =F_int(r)/(α D_xx)^2b+(4π N_Sk)^2b-4π N_Sk/r.

Since F_int(r) < 0, a skyrmion moves to the nanodisk's center along a  counter-clockwise spiral trajectory (ṙ<0, ϕ̇>0).  When it approaches the center, the angular velocity decreases since F_int(r) drops faster to 0 than 1/r. The transverse component of the spiral motion is a consequence of the skyrmion's topological charge. Like any quasiparticle with nontrivial real-space topology a skyrmion does not move parallel to the propelling force — let it be due to a current or the skyrmion–edge interaction <cit.> — and thus exhibits a skyrmion Hall effect. The corresponding transverse force is similar to the Magnus force for classical particles <cit.>.



 §.§ Formation of biskyrmions

The formation of a biskyrmion starts with two Bloch skyrmions with opposite helicities of γ = ±π/2, as stabilized in Fig. <ref>. Since both skyrmions have identical dissipative tensors D and topological charges N_Sk=+1, the equations of motion (<ref>) hold for both quasiparticles. They move on a spiral trajectory towards the nanodisk's center. However, the skyrmion-skyrmion interaction has to be considered in F⃗_int. This force changes sign when both skyrmions come close to each other, corresponding to the energy minimum in Fig. <ref>a; it is zero at a finite distance of the two skyrmions.



This behavior is well visible in the results of micromagnetic simulations (Fig. <ref>). After the writing period (relaxation after ≈ 1 ns as presented in Fig. <ref>), both skyrmions start to exhibit a spiral motion to the center as in the case of a single skyrmion. Once both skyrmions start to overlap the spiral motion stops and a biskyrmion is formed. During the formation process the topological charge remains N_Sk=+2 within an error of less than [1.5]% (the topological charge is not integer due to the discreteness of the underlying lattice).

The biskyrmion configuration is remarkably stable. The energy drops significantly once the two skyrmions are close to each other (Fig. <ref>d). This drop corresponds to the energy minimum established in the analytic model (Fig. <ref>a). In the final phase of the propagation the biskyrmion aligns with one of the crystallographic directions α = 45^∘, 135^∘, 225^∘ or 315^∘, visible as energy minimum in Fig. <ref>c (here: α = 135^∘). Once the system is close to the energy minimum F⃗_int is far from being rotational symmetric and the equations of motion derived from the Thiele equation do not hold anymore. By comparing the energy scales of Figs. <ref>a and <ref>c it becomes apparent that the reorientation of the biskyrmion plays a role after the optimal distance of the two skyrmions has been reached.
The distinguished axis of the biskyrmion can in principle be imposed by the initial positions of the skyrmions. If the two writing disks are closer together the biskyrmion can end for example at α = 45^∘ (Fig. S1 of the Supplementary information).

Starting with two skyrmions with identical helicity (i. e. by switching one of the two writing vortices) a mostly static behavior is observed (Fig. S2 of the Supplementary information). This is readily explained by the stronger repulsion of skyrmions with the same helicity in comparison to those with opposite helicities. Due to the disk geometry both skyrmions experience a force towards the center, resulting in a zero of F_int at larger distances compared to the biskyrmion system. The (final) steady state — two skyrmions without considerable overlap —  has a sizably higher energy than the biskyrmion state (cf. Figs. S2d and Fig. <ref>d). Starting from this particular configuration a biskyrmion does not form since one skyrmion would have to reverse its helicity. That, however, requires to overcome an energy barrier similar to that in Fig. <ref>b, corresponding to a Néel skyrmion state.

In Refs. <cit.> it was shown that asymmetric skyrmion tubes (e. g. in the cone phase of chiral magnetic materials) can exhibit an attractive interaction despite their equal helicities. The resulting non-overlapping pair of skyrmions is stabilized by a different mechanism and at a larger inter-skyrmion distance compared to the type of biskyrmion discussed in the present paper. Nevertheless, the skyrmion-skyrmion pair was labeled `biskyrmion' in Ref. <cit.> due to the topological equivalence (N_Sk=± 2) of both textures. Yet, both are geometrically distinct objects.



§ DISCUSSION

Here, we have demonstrated that the dipole-dipole interaction in centrosymmetric magnetic systems can stabilize isolated biskyrmions, in the sense of a pair of partially overlapping skyrmions with reversed in-plane magnetizations. We propose methods of writing two Bloch skyrmions with opposite helicities at opposite sides of a nanodisk. After spiral propagation of each skyrmion towards the nanodisk's center, an isolated biskyrmion is formed due to the attractive interaction that originates in the dipole-dipole interaction. These ingredients provide a step towards utilizing biskyrmions as carriers of information in spintronics devices. 

To the best of our knowledge, only biskyrmion lattices have been found experimentally so far <cit.>, and even for them the interpretation of the presented Lorentz transmission electron microscopy images is under serious debate. In two recent publications it was shown that tubes of topologically trivial bubbles (called type-II or hard bubbles) can appear as biskyrmion-like features when observed under an angle <cit.>. Since thin magnetic disks are considered in the present work, potential real-space images can hardly be misinterpreted, and therefore our prediction may be decisive for showing that biskyrmions can exist at all.

Besides their potential for applications, biskyrmions are also worth being investigated from a fundamental point of view. Their nontrivial real-space topology, manifested in the topological charge of N_Sk = ± 2, imposes emergent electrodynamic effects: just like conventional skyrmions, biskyrmions exhibit a skyrmion Hall effect <cit.> and a topological Hall effect <cit.>. In other words, under the effect of a spin-polarized current the biskyrmion and the current electrons are deflected into transverse directions. Our study motivates future in-depth analyses of similarities and differences of these Hall effects for skyrmions and biskyrmions.

One fundamental aspect that distinguishes biskyrmions from skyrmions is their missing rotational symmetry. The orientation of a biskyrmion could be exploited to store multiple bits per biskyrmion. Here, the underlying cubic lattice yields four energetically degenerate orientations α and allows for quaternary instead of binary logic. It is conceivable to store data in an array of nanodisks instead in a racetrack device.



§ METHODS

For the micromagnetic simulations the Landau-Lifshitz-Gilbert equation <cit.>

    ṁ⃗̇_i=-γ_em⃗_i×B⃗_i,eff+α_gm⃗_i×ṁ⃗̇_̇i̇+γ_e ϵβ[(m⃗_i×s⃗_i)×m⃗_i]

is solved using Mumax3 <cit.>.

The first term describes the precession of each magnetic moment around its instantaneous  effective magnetic field B⃗_eff^i = -δ_m⃗_i F /M_s. This field, derived from the free energy density F, takes into account exchange interaction, uniaxial anisotropy, Zeeman energy, and the demagnetization field accounting for the dipole-dipole interaction. 

The second term, introduced phenomenologically, considers that the magnetic moments tend to align along their effective magnetic field. The strength of this damping is quantified by the Gilbert damping constant α_g. 

The third term accounts for the coupling of the magnetic moments with injected spin-polarized currents. The in-plane torque coefficient ϵβ = ħ P j / 2ed_z M_s depends on the current's spin polarization P and the current density j. The spatial distribution of the injected spins {s⃗_i} is given by the vortex textures in the writing disks. An  out-of-plane torque has been neglected since it is usually small and merely rescales the external magnetic field. γ_e = 1.760 × 10^11T^-1s^-1 is the gyromagnetic ratio of an electron.

We use parameters that stabilize stray-field skyrmions, similar to those in Ref. <cit.>: exchange A = 15 pJ/m, uniaxial anisotropy in z direction K_z = 1.2 MJ / m^3, external field B⃗ = -40 mT e⃗_z (along the ferromagnetic orientation), saturation magnetization M_s = 1.4 MA / m, and Gilbert damping α_g = 0.3. The disk has a radius of 150 nm and a thickness of d_z=3 nm. For the simulations we used discrete magnetization cells of size 1 nm× 1 nm× 1 nm. The DMI has been set to zero. We checked that a small DMI constant (interfacial DMI D = 0.1 mJ/m^2) does not qualitatively affect the results of the simulations shown in Fig. <ref>; for larger values of D one of the two skyrmions would annihilate and a biskyrmion would not form.

As stated in the main text, a spin-polarized current of j⃗_writeP = 7× 10^8 A/cm^-2e⃗_z with a spin polarization orientation {s⃗_i}, given by an in-plane vortex texture

    s⃗_i=±(-y,x∓Δ x/2,0)/√((x_i±Δ x/2)^2+y^2),

writes well-defined skyrmions into the nanodisk, if the current is applied for 200 ps. The two signs correspond to the left and right writing disks, respectively. The writing disks have a radius of 25 nm each and their centers are at a distance of 120 nm. The writing disks and the process of spin polarization have not been simulated explicitly but are condensed into the spin-polarization parameter P which scales the necessary applied current density.





§ ACKNOWLEDGEMENTS

B.G. is thankful to Alexander F. Schäffer regarding discussions about Mumax and possible alternative writing approaches. This work is supported by Priority Program SPP 1666 and SFB 762 of Deutsche Forschungsgemeinschaft (DFG).



§ AUTHOR CONTRIBUTIONS STATEMENT

B.G. initiated the project and performed the calculations. B.G. wrote the manuscript with the help of J.H. All authors discussed the results and commented on the manuscript.



§ ADDITIONAL INFORMATION



 §.§ Supplementary information

Supplementary information accompanies this paper at [URL].



 §.§ Accession codes

The code Mumax3 is accessible at https://github.com/mumax/3/releases/tag/v3.9.3.



 §.§ Competing interests

The authors declare no competing interests.

