Edge instabilities and skyrmion creation in magnetic layers

We study both analytically and numerically the edge of two-dimensional ferromagnets with Dzyaloshinskii-Moriya (DM) interactions, considering both chiral magnets and magnets with interface-induced DM interactions. We show that in the field-polarized ferromagnetic phase magnon states exist which are bound to the edge, and we calculate their spectra within a continuum field theory. Upon lowering an external magnetic field, these bound magnons condense at a finite momentum and the edge becomes locally unstable. Micromagnetic simulations demonstrate that this edge instability triggers the creation of a helical phase which penetrates the field-polarized state within the bulk. A subsequent increase of the magnetic field allows to create skyrmions close to the edge in a controlled manner.


I. INTRODUCTION
The presence of the Dzyaloshinskii-Moriya (DM) interaction in ferromagnets favours a spatial twist of the magnetization leading to modulated magnetic textures like helices and skyrmion lattices, i.e., closely packed arrangements of single skyrmions. Skyrmions are spatial configurations of the magnetization with a finite topological winding number. They are not only found in the thermodynamic ground state but single skyrmions also arise as topologically protected excitations of the field-polarized state. Skyrmion couple very efficiently to spin currents so that ultralow current densities are already sufficient to drive skyrmion configurations 1,2 . Both properties, their topological protection as well as their mobility, identify skyrmions as promising candidates for elementary units of future spintronic devices, see the reviews of Refs. 3-5. A prerequisite for a skyrmion technology 6 is, however, the stabilization of magnetic skyrmion configurations at ambient conditions as well as their controlled writing and deleting. The existence of skyrmions far below room temperature has been experimentally demonstrated some years ago in the chiral magnets with space group P2 1 3 like MnSi 7,8 , FeGe 9 or Cu 2 OSeO 3 10,11 as well as in magnetic mono-and bilayers 12,13 . Recently, various groups have reported the observation of skyrmion configurations at ambient temperatures. Co-Zn-Mn alloys with the chiral space group P4 1 32 were shown to possess the typical phase diagram of other chiral magnets but with a skyrmion lattice phase stabilized at room temperature 14 . In addition, certain magnetic multilayers comprising magnetic Co and Fe atoms also show stable skyrmion configurations at room temperatures [15][16][17] . In all these systems, the DM interaction arises due to a lack of inversion symmetry. The atomic crystal of chiral magnets explicitly breaks inversion symmetry so that the DM interaction is even present in bulk samples. In magnetic multilayers, on the other hand, the DM interaction is induced by the interface between a magnetic layer and a substrate containing heavy elements.
As skyrmion configurations are topologically pro- Chain of skyrmions at the edge of a twodimensional chiral magnet that were created with the help of the magnetic field protocol shown in the inset exploiting the local edge instability of the field-polarized state. These results were obtained with micromagnetic simulations of the Landau-Lifshitz-Gilbert equation. The color code reflects the z-component of the magnetization. tected, their creation or annihilation usually requires the climb of a topological energy barrier in some way or another. It has been theoretically discussed that magnetic skyrmions can be nucleated by the local injection of a spin-polarized current 18 , local heating 19,20 , local application of a magnetic 20 or electric field 21 , and by thermal fluctuations 22 . Experimentally, it has been already demonstrated that skyrmions can be created in magnetic multilayers with the help of local currents from an STM tip 23 . Magnetic skyrmions also materialize when driving stripe domains through a geometric constriction with the help of a spin-current 16,24 . It was also shown that skyrmion crystals can be unwinded with the help of Bloch points 25 . Such Bloch points correspond to magnetic monopoles in the emergent electrodynamics 2,26 that arises when itinerant electrons locally adapt their magnetic moment to the skyrmion configurations 3,5 . Climbing the topological energy barrier and the concomitant singular Bloch point configuration of the magnetization can be avoided however by feeding in skyrmions from the edge of the sample. In fact, recent experiments on FeGe nanostripes have beautifully demonstrated already the edge-mediated nucleation of skyrmions 27 . In the present work we provide a theoretical explanation of this phenomenon. In addition, we demonstrate that with the help of a certain magnetic field protocol skyrmions can be created in a controlled manner at the edge of a magnetic monolayer by exploiting a local edge instability of the field-polarized state, see Fig. 1.
For this purpose, we investigate the edge of a single magnetic layer with DM interaction that is polarized by a perpendicular magnetic field. The boundary conditions arising from the DM interaction result in a twist of the magnetization close to the edge which was discussed independently by Wilson et al. 28 and Rohart and Thiaville et al. 29 . This surface twist was experimentally investigated in Refs. 27 and 30, and it also plays an important role in the stabilization of skyrmion configurations in thin films of finite thickness [31][32][33] . Here, we focus on the magnon excitations of the field-polarized layer, and we demonstrate that the surface twist result in spin-wave excitations that are bound to the edge of the magnetic layer. Within a continuum field-theory, we determine analytically the effective Hamiltonian of these bound magnon modes and evaluate their spectrum for various values of the effective magnetic field and the magnetic anisotropy. For vanishing H = 0, we reproduce the results of a recent numerical study of these edge excitations 34 . Moreover, for finite H and intermediate values of the magnetic anisotropy the edge magnons become locally unstable and condense at a finite momentum. It is this local instability that facilitates the skyrmion creation.
The outline of the paper is as follows. In sec. II we determine the magnetization at the edge of a single layer, we derive the magnon Hamiltonian and discuss its eigenmodes that are bound to the edge of the magnetic layer. We summerize the instabilities of the field-polarized state in sec. II E with a particular focus on the local instability triggered by the bound states that become unstable at a finite transversal momentum. In sec. III we demonstrate with the help of micromagnetic simulations that this latter instability triggers the formation of a helical phase, and, in addition, how this instability can be used to create skyrmions. We conclude in sec. IV with a short discussion.

II. MAGNONS AT THE EDGE OF A CHIRAL MAGNET
A. The Free energy functional For our study, we consider the one-dimensional edge of a two-dimensional magnetic layer. In the following, when we speak about the bulk of the system we refer in fact to the interior of this two-dimensional layer. For sufficiently low temperatures, one can neglect (thermal) amplitude fluctuations of the spins and the local direction of the magnetization can be described by the unit vectorn. All modes discussed in the following including the boundary states at the edge are characterized by length scales much larger than the lattice spacing, so that a description in terms of a continuum model is applicable.
The free energy functional of the magnetic layer, F = d 2 rF, reads with α = x, y and i = 1, 2, 3, A is the stiffness, and K is the magnetic anisotropy with an easy-axis and easyplane anisotropy for K > 0 and K < 0, respectively. Note that for a two-dimensional layer the dipolar interactions are effectively accounted for by a renormalized anisotropy K 35 . The third term is the Zeemann energy with an external field H applied perpendicular to the plane, H = Hẑ with H > 0, and the magnitude of the local magnetization M (integrated over the z direction). There a two different types of Dzyaloshinskii-Moriya (DM) interactions F DM . For chiral magnets like MnSi or FeGe where the inversion symmetry is broken by the atomic crystal the DM interaction reads 36 where iαj is the Levi-Civita symbol, while for interfaceinduced DM interactions relevant, e.g., for Co films grown on heavy-element layers (e.g. Pt or Ir), we have 37,38 with α = x, y and i, j = 1, 2, 3. It has been pointed out in Ref. 39 that for two-dimensional systems both types of DM interactions are actually equivalent because they can be transformed into each other by rotatingn by π/2 around the z-axis, i.e., byn x →n y ,n y → −n x . We will use in the main part of this work the formulae appropriate for F chiral DM and assume D > 0. Our results are straightforwardly generalized to interface-induced DM interactions with the help of the above symmetry.
The dynamics of the magnetization is governed by the equation of motion with the gyromagnetic ratio γ = gµ B / > 0 and the effective magnetic field obtained by the functional derivative B eff = − 1 M δF/δn where F = d 2 rF. For our micromagnetic simulations, we furthermore add Gilbert damping to these equations, see Ref. 40. The typical momentum, time and energy scales are given by In the main part of the paper, we will measure length, time and energy in dimensionless units by denoting the variables with a tilde, e.g.
Moreover, it is convenient to introduce the dimensionless magnetic field h and the dimensionless anisotropy κ, As we are interested in the properties of the edge of the layered sample, we have to supplement the equations of motion (4) by a boundary condition. We assume the edge parallel to the x-axis. By varying the free energy F = ∞ −∞ dx ∞ 0 dyF, one obtains for chiral magnets described by Eq. (2) the boundary condition 30 ∂ yn − Qŷ ×n = 0 at y = 0 (8) while for interface-DM interactions (3) one finds the boundary condition 29 ∂ yn + Qx ×n = 0 at y = 0.
A subtle question is whether the continuum approach advocated above is valid for real materials. In general, the chemistry and therefore the exchange couplings, the anisotropies, and DM interactions at the edge will differ from the ones within the layer. These changes do, however, affect our results only weakly if (i) the modification of the chemical properties is limited to distances of a few lattice constants a away from the edge and (ii) spin-orbit coupling effects are weak. The latter condition implies that the typical length scale on which the magnetization varies is much larger than the lattice constant a and, furthermore, that anisotropy terms remain small compared to exchange interactions even at the surface. Technically, one can check for the importance of surface modification by adding surface, i.e., edge terms to the continuum Lagrange density, e.g., of the form K s δ(y)n 2 3 , A s δ(y) (∂ yni ) 2 , or D s δ(y)n 1 ∂ xn2 with coupling constants K s , A s , D s . The δ-function takes into account that only the edge is affected. The importance of such terms can be evaluated by doing power-counting in the strength of spin orbit coupling λ SO where we assume such that the dimensionless parameters h and κ are of order 1. Under this condition all length scales are proportional to 1/λ SO . As δ(y) = λ SO δ(yλ SO ) all extra terms discussed above are suppressed by the small factor λ SO compared to the bulk contributions. This argument does nominally not hold for the surface term δ(y)n∂ yn which by powercounting is as important as the bulk contributions but this term is identical to zero as ∂ yn 2 = 0. For a discussion of surface terms for a three-dimensional magnet see Meynell et al. 30 .

B. The magnetization at the edge of the layer
We first consider the static magnetization close to the edge assuming that the magnetic layer is in a polarized state withn =Ĥ =ẑ within the bulk of the layer. Close to the edge of the two-dimensional chiral magnet, the magnetization has to twist due to the boundary conditions (8) or (9). An analytic solution for this twist has already been derived by Meynell et al. 30 for a system without magnetic anisotropy K = 0. Below, we generalize their result including a finite K. In order to simplify notations we use the rescaled coordinates of Eq. (6).
The translational invariance along the edge ensures thatn is only a function ofỹ, i.e., the distance from the edge. We consider a chiral magnet described by Eqs. (2) and (8) where the magnetization twists perpendicular to the propagation direction similar to a Bloch-type domain wall. Hence we choosê as an ansatz for the magnetization. The equation of motion (4) in the static limit and the boundary condition of Eq. (8), respectively, then reduce to Moreover, θ(ỹ), θ (ỹ) → 0 forỹ → ∞ as we assume the magnetizationn =ẑ to be aligned with the field within the bulk of the layer.

Kink soliton of the double sine-Gordon model
The equation of motion Eq. (11) is just the Euler-Lagrange equation of the double sine-Gordon model 41,42 . It turns out that the magnetization at the edge can be obtained by placing a kink soliton of this model close to the edge. The analytical expression for this kink soliton is known in closed form 41,42 and we shortly discuss its derivation in the following. As the derivative at the edge is positive, θ (0) > 0, we will look for a kink that interpolates θ from −2π to 0 for increasingỹ. The first integral of the equation of motion is obtained in a standard fashion by multiplying Eq. (11) with θ , integrating and using the boundary conditions within the bulk of the layer, θ(∞) = θ (∞) = 0, The kink is obtained by solving Eq. (13) for θ (ỹ) > 0. Integrating Eq. (13) by separation of variables yields the kink soliton whereỹ 0 is an integration constant that specifies the position of the kink. For later convenience, we will also discuss the energy of the kink when it is placed deeply within the bulk,ỹ 0 → ∞. The kink energy per length L x is obtained by integrating where F kink is the free energy density evaluated for the kink, and we have subtracted the energy density of the field-polarized state, F FP , There exists a critical magnetic field h cr kink for which the kink energy vanishes, This critical field is defined within the range −1 < κ < π 2 /4 where it monotonically decreases from h cr kink = 1 to zero. For κ = 0 one recovers h cr kink = π 2 /16 (Ref. 38). Whereas for larger fields the kink is an excitation of the field-polarized state, the commensurate-incommensurate transition can take place at h cr kink , and kinks become energetically preferred resulting in a soliton lattice. However, as the kink is a topological excitation it does not correspond to a local but rather a global instability of the field-polarized state for h ≤ h cr kink .

Edge magnetization
The full magnetization profile is now easily obtained by placing the kink Eq. (14) close to the edge so that the boundary condition, θ (0) = 1, is fulfilled. This is achieved with the integration constant The center of the kink is positioned outside the sampleỹ 0 < 0. Note that there is also a solution of the boundary condition with a positiveỹ 0 but it yields a state with larger energy. The result of Eq. (14) together with Eq. (17) determines the magnetization close to the edge of the field-polarized magnetic layer.
In the limit of vanishing anisotropy, κ = 0, we recover the results of Ref. 30, We checked the result for various combinations of h and κ also numerically by simulating the damped Landau-Lifshitz-Gilbert equation and find that numerics are in excellent agreement with the analytic expression. In Fig. 2 we show the analytic result for the z-component of the magnetization in the proximity of the edge.

C. Effective magnon Hamiltonian
Having established the profile of the magnetization, we now consider the fluctuations around the saddle-point solution defined by Eqs. (10), (14) and (17) following the approach of Ref. 43. We introduce a local basiŝ with the angle θ = θ(ỹ) given by Eq. (14). We use the following parametrization of the fieldn, withê ± = 1 √ 2 (ê 1 ±iê 2 ). A change of phase of the complex magnon wavefunction ψ = ψ(x,ỹ,t) naturally encodes the precession of the magnetization.
Expanding the Landau-Lifshitz equation (4) up to linear order in the wavefunction ψ and ψ * , we obtain the wave equation for the two-component spinor ψ = ψ ψ * , τ i are Paulimatrices, and H is an effective Bogoliubov-deGennes Hamiltonian.
It can be split into the bulk contribution H 0 and the edge potential V = V (ỹ) that vanishes V → 0 forỹ → ∞, After performing a Fourier transformation of the spinor wavefunction, ψ(q x ,ỹ,t) = dxe −iqxx ψ(x,ỹ,t), for thẽ x direction along which the problem is translationally invariant, we obtain whereq x is the dimensionless wavevector alongx. The magnon gap within the bulk of the layer is given by For ∆ b < 0, the field-polarized state within the bulk is locally unstable with respect to a tilt of the magnetization away from the field axis. The potential reads with the angle θ(ỹ) given in Eq. (14). We have used the first integral (13) to simplify the potential and, in particular, to eliminate the explicit dependence on dimensionless magnetic field h. Moreover, an explicit calculation shows that the boundary conditions for the spinor is simply given by The bosonic Bogoliubov-de Gennes Hamiltonian H, Eq. (25), possesses scattering states that extend into the bulk of the layer as well as magnon states that a bound to the edge by the potential V . In the following, we will concentrate on these magnon edge modes.

D. Bound magnon edge modes
We will look for eigenstates φ = φε ,q (ỹ) that are localized at the edge of the chiral magnet, which requires the dimensionless energyε to have values within the range, 0 ≤ε < ∆ b +q 2 x . The localized eigenstates obey the stationary wave equation 43 Hφ =ετ z φ, with the normalization condition From the spinor-wavefunction φ T = (φ 1 , φ 2 ) then follows the corresponding time-dependent magnon wavefunction in the parametrization of Eq. (23), that is We numerically search for bound state solutions using the shooting method. Starting from the boundary condition at the edge φ T (0) = c 1 (cos α, sin α) and φ (0) = 0 where c 1 is an arbitrary constant that will be fixed afterwards by the normalization condition (31), we vary the energyε and the parameter α until we find a solution of Eq. (30) that is bound to the edge so that it decays for y → ∞. The asymptotic decay of the localized wavefunction directly follows from H 0 of Eq. (26), with constants c 2 and c 3 . The problem (30) corresponds to an effective one-dimensional wave equation, and its bound states can be labelled by a discrete quantum number n = 0, 1, 2, ... that count the nodes of the wavefunction. Moreover, the solutions for given momenta along the edge q x define dispersing eigenenergies ε n (q x ) = ε DMεn (q x ) withq x = q x /Q for the magnon bound states that are discussed in the following.

Lowest-energy magnon edge modes
The dispersion of magnon edge modes with lowest energy ε 0 (q x ), i.e., quantum number n = 0 and momentum q x along the edge are shown for various values of the dimensionless magnetic field h and κ = 0 in Fig. 3 as well as for various values of the dimensionless magnetic anisotropy κ and h = 0 in Fig. 4. We have rescaled the energy by the dimensionless bulk gap ∆ b = h + κ and the momenta by √ ∆ b such that the bottom of the bulk continuum, ∆ b +q 2 x , is given by a single black line for all parameters. Our results are in agreement with the numerical study by Garcia-Sanchez et al. in Ref. 34.
Generally, we always find a range of momenta q x where magnon modes exist that are bound to the edge. Interestingly, the spectrum is chiral, i.e., it is not symmetric around q x = 0 as the presence of the magnetic field breaks all q x → −q x symmetries at the boundary. Especially for larger magnetic fields h the group velocity ∂ε 0 (q x )/∂q x is mostly positive, which implies that magnons travel along the edge in a preferred direction 34 , for example, after they have been excited locally by a laser pulse. Note that the modes with the reversed dispersion ε n (q x ) → ε n (−q x ) are found at opposite edges of the layer. Theoretically speaking, the dispersion could be reversed at the same edge by reversing the sign of the gyromagnetic ratio γ → −γ, e.g., by considering holes instead of particles. However, reversing the sign of the DM interaction, D → −D, modifies the edge magnetization but eventually leaves the dispersion ε n (q x ) invariant in contrast to the claim of Ref. 34. This also applies for interfacial DM interaction. While for κ = 0 (Fig. 3) bound states exist only for q x > 0, one obtains edge magnons at q x = 0 for h = 0 and κ > 1.005. In the latter case an homogeneously oscillating magnetic field can be used to excite a left-moving edge mode. In the presence of disorder at the edge, i.e., when the momentum q x parallel to the edge is not conserved, also right-moving modes can be excited by such an ac field.
With lower magnetic fields the minimum of the spectrum, is shifted down to lower energies and to larger wave vectors. At a critical dimensionless magnetic field h c = 0.4067 at κ = 0 the minimum of the spectrum goes to zero energy, ∆ e = 0, at a finite momentum q x . The system therefore experiences a local instability at the edge which we discuss in further detail in Sec. II E. Similarly, we find such a local edge instability for κ c = 1.005 and h = 0.

Higher-order magnon edge modes
For larger values of q x , the gap between the lowest bound magnon mode ε 0 (q x ) and the bulk continuum increases approximately linear. This dependence can also be estimated from the potential, Eq. (28), which grows linear in q x for large q x . As the potential becomes deeper, higher-order bound states, ε n (q x ) with n > 0, arise at the edge which are characterized by one (or more) nodes in their wave function. The spectrum of the bound states with zero and one node is shown in Fig. 5 for the magnetic field h = 0.44.

E. Phase diagram and instabilities of the metastable field-polarized state
In the following we discuss the phase diagram of the magnetic layer as a function of magnetic field h and magnetic anisotropy κ, see Fig. 6. First, we review the thermodynamically stable phases, and, afterwards, we focus on the metastable field-polarized phase and its global and local bulk and local edge instabilities.

Thermodynamic phase diagram
The bulk thermodynamic phase diagram of a magnetic layer in the presence of a perpendicular magnetic field has been studied in Refs. 39, 44-46. For a twodimensional chiral magnet, there exist three stable thermodynamic phases for values of κ that are of interest here: the field-polarized state (FP), the skyrmion crystal (SkX), and the chiral soliton lattice (CSL), i.e., a helix that in general possesses higher harmonics 47 . The black dashed lines in Fig. 6 show the phase boundaries that were given by Wilson et al. 44

Global, local bulk and local edge instabilities of the metastable field-polarized state
Importantly, the transition between the field-polarized state (FP) and the skyrmion crystal (SkX) corresponds to a global instability as the two states belong to different topological sectors. If the magnetic field is reduced adiabatically at low enough temperatures to a value within the SkX regime, there is the possibility that the magnetic layer remains in a metastable field-polarized state. As the field is reduced further, another global instability of the metastable FP state is encountered at h cr kink within the range −1 < κ < 1.9 that is shown as the dashed green line in Fig. 6. Here the energy of a kink soliton vanishes, see Eq. (16), triggering the formation of a chiral soliton lattice. The phase boundary between FP and CSL state within the range 1.9 < κ < π 2 /4 is also defined by h cr kink . The field-polarized state is stable or metastable in the whole white regime of Fig. 6. It becomes locally unstable however when one of its magnon excitations reaches zero energy upon decreasing the magnetic field h. This instability can occur either in the bulk or at the edge. For κ < −0.61, first the bulk becomes locally unstable at h = −κ (grey solid line) where the bulk gap vanishes, ∆ b = 0, see Eq. (27). For −0.61 < κ < 1.005, on the other hand, the gap of the magnon edge modes vanishes first, ∆ e = 0, see Eq. (33), at the red solid line. This edge instability occurs at a finite transversal momentum q x , see Fig. 3 and Fig. 4.
Finally, the blue solid line in Fig. 6 indicates where the bulk and edge gap have equal size. Below that line the gap ∆ e of the edge magnons is smaller than the gap ∆ b of the bulk magnons, and the spin wave excitation with lowest energies are located at the edge of the sample. Within this regime, for frequencies ∆ e ≤ ω < ∆ b only edge magnons of the stable or metastable FP state are thus excited.

III. EDGE INSTABILITY AND CREATION OF SKYRMIONS
In the following we demonstrate that the local edge instability of the metastable field-polarized state at the red line in Fig. 6 triggers the formation of a helical state (CLS) although in some region of the phase diagram the skyrmion crystal (SkX) is in fact thermodynamically more stable. Moreover, we show that this edge instability can be exploited to create skyrmions in a controlled manner. In order to investigate the evolution of the magnetic state in the regime where the field-polarized state is locally unstable, we have performed micromagnetic simulations using the Landau-Lifshitz-Gilbert (LLG) equation at T = 0, for details see Ref. 40. Typically we use α = 0.1 or 0.4 as a damping constant in our simulations. We simulated a two-dimensional stripe with open boundary conditions in one and periodic boundary conditions in the perpendicular direction.
We focus on the range of magnetic anisotropies where the edge magnons locally destabilize the FP state, i.e., −0.61 < κ < 1.005 for the two-dimensional system, see Fig. 6. The following protocol for the time-dependent magnetic field allows for a controlled creation of a skyrmion chain close to the edge The initial field value h i is located above the upper black dashed line in Fig. 6 so that we start with a magnetization which is field-polarized. At t = 0 the field is reduced h 0 < h i below the red line in Fig. 6 so that ∆ e < 0. Finally, at time t f we increase the field again to a value, h f > h cr kink , that is located above the green dashed line in Fig. 6. The result of the micromagnetic simulation for such a protocol at κ = 0 is shown in Fig. 7.
Initially, the magnetization is fully polarized except close to the lower edge of Fig. 7(a) that shows the surface twist of Eq. (14). As the field is lowered for t > 0, the edge magnon becomes soft at a finite transversal momentum q x and destabilizes the magnetization whose time evolution is shown in Fig. 7(b)-(d). To trigger the edge instability in the numerics, we explicitly broke translation symmetry by introducing a tiny perturbation by canting one spin at the right-hand side of the boundary by 1%. First, a periodic modulation of the edge spins grows in amplitude -the edge magnon with negative energy becomes macroscopically occupied at finite momentum q x . This state evolves smoothly into a helical state which penetrates into the field-polarized state of the bulk. The interface between the helical and polarized phase is thereby described by a chain of merons, i.e., half-skyrmions with winding number 1/2. As a function of time, the interface moves further and further into the field-polarized state of the bulk. For this simulation the intermediate field value h 0 of Eq. (34) was chosen to be located within the regime where the skyrmion crystal (SkX) is thermodynamically stable. Nevertheless, the local edge instability prompts the formation of a helical state which is metastable in this case.
In a second step, the magnetic field is again increased to a value h f at t = t f , and, as a consequence, the helical phase is pushed back towards the edge. For this to happen, it is required that h f > h cr kink so that the fieldpolarized state is energetically favoured compared to the CSL state and can exert a positive pressure on the interface. For the simulation in Fig. 7(e)-(g), we have chosen h f = h i . Remarkably, the initial state is however not recovered. An interesting local dynamics governs the fate of the helical fingers as they are pushed towards the edge. When the merons, i.e., the half-skrymions approach the boundary, each of them pulls a second meron out of the edge. Both combine to a skyrmion which gets repelled from the edge by the surface twist of the magnetization. As a result, one obtains a chain of equally spaced skyrmions. Note that h f = h i is here located within the regime where the FP phase is thermodynamically stable. The interaction of the receding meron with the spin configuration at the edge thus results in a final state, Fig. 7(g), that possesses a larger energy than the initial field-polarized state of Fig. 7(a). The spin configuration at the boundary with its surface twist apparently acts as a repulsive potential for the merons hindering them to leave the sample so that the field-polarized thermodynamic groundstate becomes dynamically inaccessible. The systems is thus trapped in a metastable state containing a chain of skyrmions.
The advantage of the protocol (34) is that it is extremely robust. The final state is, for example, completely independent of the time t f for which the magnetic field is lowered. The precise field values and the size of the damping term is also not important. Another important aspect of this protocol is that during the process, the magnetic configuration is always smooth and singular intermediate spin configurations, i.e., Bloch points are not needed. This is in contrast to processes where skyrmions are created within the bulk 21,25 .
The second part of our protocol, the creation of a chain of skyrmions by pushing the helical state towards the edge with the help of a magnetic field has been realized recently for FeGe nanostrips in a beautiful experiment by Du et al. in Ref. [27]. In this experiment, the system was first prepared in a helical state. After an increase of the magnetic field, helices ultimately reconstruct yielding skyrmions. For nanostrips with a small width, the skyrmions recombined in the center of the strip while for larger width precisely the same process was observed that we discuss above.
The same protocol can also be used to create a single skyrmion by reducing the magnetic field however only within a confined region close to the edge as indicated by the dashed box in Fig. 8(a). Experimentally, one can use, e.g., a magnetic tip to change locally the magnetic field. Previous micromagnetic simulation by Koshibae et al. 20 have already demonstrated that pulses of magnetic field in a small area can trigger the creation of a skyrmion close to the edge. Here, we identify the edge magnon instability as the underlying principle of this phenomenon. As the field is locally reduced below the edge magnon in- stability, a meron is protruding into the system as shown in Fig. 8(b)-(d). As the magnetic field is increased again a single skyrmion forms, see Fig. 8(e) and (f). In panel (g) we show the time evolution of the total skyrmion winding number. As the field is decreased at time t = 0, the winding number increases and saturates to a value close to 1/2 reflecting the presence of a single meron. As the field is increased again at t f = 90t DM , the winding number increases assuming a value close to one when the skyrmion formation is completed.

IV. DISCUSSION
The Dzyaloshinskii-Moriya (DM) interaction imposes boundary conditions on the magnetization 28,29 which results in a reconstruction of the magnetization profile close to surfaces. For the field-polarized state this just translates to a twist of the magnetization along or perpendicular to the surface normal depending on the type of DM interaction. We have demonstrated that this surface twist can act as an attractive potential on the spin wave excitations leading to magnon modes that are bound to the edge of the sample. The energy of these edge modes as a function of momentum transverse to the edge has been calculated and is shown in Figs. 3 and 4. These bound magnons become soft at a finite momentum thus locally destabilizing the metastable field-polarized phase, and this edge instability triggers the formation of a helical state as shown in Fig. 7. We have shown by micromagnetic simulations that the process of helix formation via the edge instability is dynamically irreversible allowing for the creation of skyrmions close to the edge.
In the present work, we confined ourselves to a twodimensional monolayer where the edge instability corresponds to a local instability of the field-polarized state. We now shortly comment on the modifications expected for a magnetic film of finite thickness. In this latter case, the field-polarized state becomes locally unstable upon decreasing the magnetic field with respect to the formation of a conical phase 44,45 . Importantly, this instability occurs before the local edge instability discussed in this paper so that our results cannot be directly applied to thin films. However, the conical phase should also exhibit its own edge instability whose analytical description is however more involved due to the coupling of the helical modulation and the surface twist. Nevertheless, we expect the edge instability of the conical phase to possess a similar character also allowing for the creation of skyrmions as discussed above. We thus believe that our results basically explain the mechanism of the edge-mediated nucleation of skyrmions experimentally observed by Du et al. 27 in thin films of FeGe.
Whereas we focused here on the properties of the field-polarized state, the surface reconstruction is also expected 31-33 for the other thermodynamically stable phases, i.e., the skyrmion crystal and the helix. It is an interesting open question how it influences the spin wave spectrum and whether or not bound magnon edge modes also exist in these other phases. At least for the skyrmion crystal phase, magnon edge modes are indeed expected but for a very different reason. Magnons experience an emergent orbital magnetic field when they scatter off a topological skyrmion configuration 3,43,48 . In a magnetic skyrmion crystal this should give rise to a topological magnon band structure characterized by finite Chern numbers and the concomitant topological edge modes 49 .